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*8  SUPPLEMENTARY  NOTES 


1**  KEY  WOROS  (Continue  on  reveree  aide  it  neceeeety  and  identity  by  btocl t  number) 

Unsteady  Transonic  Flow  Finite  Elements  Method  of  Weighted  Residuals 

Steady  Transonic  Flow  Potential  Flow  Least  Squares 


*°  Y  RACT  (Continue  on  reverse  aide  it  net:  eeaery  and  identity  by  block  number )  I 

A  finite  element  program  is  described  for  computing  steady  and  unsteady 
(oscillatory  and  transient)  transonic  flows  over  thin  airfoils  by  solving  directly 
the  unsteady,  nonlinear  transonic  potential  equation  based  on  small  disturbance 
theory.  The  present  numerical  algorithm  is  developed  using  the  concept  of 
finite  elements  in  conjunction  with  the  least  squares  method  of  weighted  resi¬ 
duals  applied  to  both  space  and  time.  The  basic  element  presently  used  is  a 
product  of  an  element  in  space  and  an  element  in  time.  The  former  has  a  cubic 


DD  ,  " 


EDITION  OF  I  NOV  «S  IS  OBSOLETE 


o^r 


UNCLASSIFIED _ _ 

SECURITY  CLASSIFICATION  OF  THIS  PAOf 


cZ/C/00' 


UN  CLASS]  !•  i  FI) 


>tcuwnv  ci  ^sihcahon  of  this  pageix/imi  (>•<• 


■xpansion  inside  each  element,  while  the  latter  is  a  quadratic  Lagrangian  ele¬ 
ment.  Fur  each  time  step,  tlie  finite  element  discretization  in  both  space  and 
time  results  in  a  recurrence  relationship  in  the  form  of  a  banded  system  of 
algebrnu  equations,  whirl  is  solved  by  Gaussian  elimination.  The  embedded 
shucks  arc  smeared  and  a  matching  scheme  for  computing  effectively  flow  over 
lilting  airfoils  is  also  incorporated  in  the  program.  The  preSent  computer 
program  is  composed  of  two  parts;  the  first  part  (designated  as  UTRANL-I) 
generates,  from  a  limited  number  of  input  cards,  the  necessary  mesh  informa¬ 
tion  and,  if  desired,  produces  a  CALCOMP  mesh  plot;  the  second  part 
1  KAVl  -11)  cairies  out  the  analysis  and  displays  the  pressure  coefficients 
along  the  chordline  on  printer  plots.  Two  sample  cases  of  flow  over  a  NACA 
o4A  410  and  a  NACA  64A  006  airfoils  are  given  to  demonstrate  the  applicability 
and  usage  of  the  program.  The  solution  procedures  are  found  to  be  quite  effl¬ 
uent  and  accurate,  permitting  the  aerodynamic  forces  to  be  calculated  to  engi¬ 
neering  accuracy  in  less  than  ten  minutes  CPU  time  on  a  CDC  6600  computer 
for  the  most  time  consuming  case  among  all  those  studied. 
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This  report  was  prepared  by  personnel  in  the  Engineering  Sciences 
Section  of  the  Lockheed  Missiles  &  Space  Company,  Inc.,  Huntsville  Research 
&  Engineering  Center,  Huntsville,  Alabama,  for  the  Air  Force  Flight  Dynamics 
Laboratory,  Wright-Patterson  Air  Force  Base,  Ohio.  The  computer  programs 
were  developed  under  Project  1370,  "Dynamic  Problems  in  Flight  Vehicles," 
Task  137004,  "Design  Analysis,"  Contract  F33615-75-C-3  125.  Capt.  Gerald 
Van  Keuren,  AFFDL/FBR,  is  the  Air  Force  Project  Engineer. 

S 

S.T.K.  Chan  was  the  principal  investigator  for  the  study,  and  H.C.  Chen, 
as  co -investigator ,  contributed  to  the  development  of  the  numerical  solution 
method  and  developed  most  of  the  necessary  computer  programs.  The  study 
was  conducted  under  the  supervision  of  M.  R.  Brashears  and  B.H.  Shirley. 

The  authors  submitted  this  report  in  March  1976  as  an  AFFDL  technical 
report  to  cover  research  performed  from  July  1975  to  March  1976.  The  theory 
and  numerical  solution  method  for  the  present  computer  program  is  documented 
as  AFFDL-TR-76-49,  Vol.I. 
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SECTION  I 
INTRODUCTION 


In  recent  years,  significant  progress  has  been  made  toward  developing 
a  useful  method  for  predicting  steady  transonic  airloads  over  airfoils  and  to 
some  degree  for  finite  wings.  Despite  this  progress,  however,  very  few  satis¬ 
factory  methods  have  evolved  for  calculating  both  steady  and  unsteady  transonic 
flows,  as  evidenced  in  a  survey  paper  by  Bland  (Ref  1).  Only  recently  have 
numerical  solutions  via  the  finite  difference  technique  been  presented  for 
two-dimensional  airfoils  executing  harmonic  motion.  More  recently,  the 
finite  element  method  has  also  been  applied  successfully  to  solve  transonic 
flow  problems  (Refs.  2  and  3). 

In  the  present  study,  the  finite  element  technique  is  extended  to  compute 
both  steady  and  unsteady  transonic  flows  over  lifting  airfoils,  based  on  small 
perturbation  theory.  Unlike  most  existing  techniques  for  solving  the  small 
disturbance  potential  equation,  the  present  approach  solves  directly  the  un¬ 
steady,  nonlinear  transonic  flow  equation,  with  both  time  derivative  terms 
retained.  Thus  the  present  algorithm  can  be  applied  to  compute  a  much 
wider  class  of  transonic  flow  problems,  including  steady,  oscillatory  or 
transient  solutions,  either  with  or  without  angle  of  attack.  For  oscillatory 
flow,  no  assumption  is  made  regarding  the  oscillating  frequency,  nor  is  the 
unsteady  perturbation  necessarily  small  compared  to  the  mean  steady 
solution. 

The  present  numerical  algorithm  was  developed  using  the  concept  of 
finite  elements  in  conjunction  with  the  method  of  weighted  residuals.  With 
the  present  approach,  the  embedded  shocks  are  smeared.  Also,  a  patching 
technique  has  been  developed  to  match  the  finite  element  solution  constructed 
in  a  moderately  large  domain  with  an  asymptotic  solution  for  the  far  field  to 
avoid  the  necessity  of  using  a  very  large  domain  otherwise  needed  in  com¬ 
putations.  The  basic  element  presently  used  is  the  product  of  an  element 
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m  spat  e  and  an  element  in  time.  The  former  has  a  cubic  expansion  inside 
the  element,  with  nodal  unknowns  representing  the  perturbed  velocity  poten¬ 
tial  and  the  two  perturbed  velocity  components;  while  the  latter  is  a  quadratic 
Lag rangian  element.  The  resulting  system  of  algebraic  equations,  which  relates 
the  current  time  step  solution  to  solutions  at  two  previous  time  steps,  is  banded 
and  can  be  solved  conveniently  by  Gaussian  elimination. 


Several  aspects  of  the  present  numerical  algorithm  are  discussed  here¬ 
in.  The  small  disturbance  potential  equation  with  associated  boundary  condi¬ 
tions  and  the  related  secondary  unknowns  are  summarized  in  Section  II.  The 
numerical  procedures  employed  in  the  present  approach  are  discussed  in 
Section  III.  More  detailed  discussion  on  the  theory  and  the  numerical  solution 
method  is  presented  in  the  first  volume  of  the  present  report.  In  Section  IV, 
the  two  parts  of  the  computer  code,  namely,  UTRANL-I  and  UTRANL-II  are 
described  separately,  in  the  aspects  of  scope  and  flow  chart,  description  of 
variables,  subroutines  used,  and  finally  input  and  output.  Two  sample  cases 
are  given  in  Section  V  to  demonstrate  how  to  use  these  computer  programs, 
which  are  listed  in  the  Appendixes. 


SECTION  II 

THEORY  -  ASSUMPTIONS  AND  BASIC  EQUATIONS 


The  objective  of  the  present  study  was  to  develop  an  efficient  and  accu¬ 
rate  numerical  algorithm  for  the  analysis  of  steady  and  unsteady  (oscillatory 
and  transient)  transonic  flows  over  thin  airfoils.  With  this  objective  in  mind, 
che  formulation  was  therefore  based  on  the  small  disturbance  but  nonlinear 
transonic  potential  equation  for  inviscid,  compressible  fluid.  The  embedded 
shock  is  assumed  to  be  weak  and  boundary  layer  effects  are  neglected.  With 
these  assumptions,  the  transonic  flow  problems  under  consideration  can  be 
stated  in  the  following  mathematical  form. 


Differential  Equation 


2M2  $  ,  -  M2 

oo  ,  xt  oo 


(1) 


Boundary  Conditions 

•  Vanishing  of  disturbance  at  the  far  field,  that  is,  at  infinity 


V<£  =  0 


(2) 


•  Flow  tangency  condition  on  the  airfoil  surface, 


•  Unsteady  Kutta  condition  in  the  wake  to  ensure  pressure  being 
continuous  in  the  vortex  sheet,  namely, 


(4a) 


and 


(4b) 


In  the  above,  <j>  =  perturbed  velocity  potential  function,  =  freestream  Mach 
number,  Y  =  ratio  of  specific  heats,  taken  to  be  1.4  for  air,  g  =  geometry  of 
the  airfoil  with  angle  of  attack  included,  6  =  amplitude  of  oscillation  and  h  = 
function  describing  the  airfoil  oscillation.  Equation  (1)  is  in  dimensionless 
form  and  the  x-axis  is  aligned  with  the  undisturbed  flow  direction.  The  dimen 
sional  (with  primes)  and  nondimen sional  quantities  are  related  by 


x 


U 

t'  -2-  and  <f> 
c 


4' 

U  c 
00 


(5) 


where  c  and  U  are  the  characteristic  length  and  speed,  which  are  currently 
00 

taken  as  the  chord  length  and  the  freestream  speed,  respectively.  Conse- 

Q 

quently,  the  characteristic  time  is  — — ,  which  is  the  time  needed  for  flow  at 

oo 

freestream  speed  to  travel  one  chord  length. 

In  actual  computations,  because  only  a  finite  domain  can  be  considered 
in  the  analysis,  an  appropriate  asymptotic  solution  instead  of  Eq.  (Z)  must  be 
imposed  on  the  outer  boundary  of  the  finite  element  mesh.  The  asymptotic 
solution  for  a  two-dimensional  airfoil  in  steady  flow  has  been  derived  by 
Klunker  (Ref.  4)  as 

*  =  £  [z  sen(y) ttan"1 
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The  first  term  on  the  right-hand  side  corresponds  to  a  free  vortex  and  repre¬ 
sents  the  lifting  effects;  the  second  term  corresponds  to  a  source  distribution 
whose  strength  is  related  to  the  airfoil  thickness  distribution  g(x);  and  the 

third  term,  which  arises  from  the  nonlinearity  of  the  flow  equations,  has  the 

2 

form  of  a  doublet  with  its  strength  given  by  the  local  value  of  u  and  is  to  be 
integrated  over  the  flow  domain  under  consideration.  Klunker  also  showed 
that  the  contributions  from  the  thickness  and  the  area  integrals  are  of  higher 
order  effects.  For  these  reasons,  only  the  first  term  is  used  currently  for 
the  far  field  in  computing  flow  over  lifting  airfoils. 

Figure  1  shows  the  flow  field,  together  with  the  unsteady  transonic 
governing  equation  and  corresponding  boundary  conditions,  in  a  finite  domain 
used  in  computations.  Numerical  experimentations  indicate  that  a  flow  region 
with  H  >2.5c  is  generally  required  to  yield  solution  with  adequate  accuracy. 
For  flow  with  higher  freestream  Mach  number,  the  flow  region  should  be  ex¬ 
tended  somewhat  to  account  for  effects  from  the  enlarging  supersonic  flow 
pockets . 


Once  the  flowfield  solution  in  terms  of  the  perturbed  velocity  potential 
has  been  obtained,  all  secondary  unknowns  can  be  calculated  subsequently. 


These  include 


(U2  -  V2  -  20  U  )  + 
2  '  00  ,  t  00 


(8) 
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Asymptotic  Solution 


Figure  1  -  Schematic  of  Flow  Field 
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(10) 


In  the  above,  =  1 ,  the  normalized  freestream  speed, 
sound,  M  =  local  Mach  number,  V  =  the  total  velocity,  p 
p  =  stagnation  pressure  and  Cp  =  pressure  coefficient. 


a  =  local  speed  of 
=  local  static  pressure, 


SECTION  III 

NUMERICAL  SOLUTION  PROCEDURES 


The  concept  of  finite  elements  in  conjunction  with  the  least  squares 
method  of  weighted  residuals  (MWR)  is  the  basis  of  the  present  numerical 
procedures  in  solving  the  small  disturbance  unsteady  transonic  flow  equation. 
The  use  of  the  least  squares  formulation  ensures  that  the  resulting  matrix 
is  positive  definite  and  well  conditioned.  The  numerical  procedures  including 
the  finite  element  formulation,  element  description,  the  imposition  of  boundary 
conditions,  and  time  marching  procedures  are  described  in  the  following  sub¬ 
sections  . 

1.  FINITE  ELEMENT  FORMULATION 

The  finite  element  method,  in  conjunction  with  the  least  squares  method 
of  weighted  residuals,  is  used  herein  to  solve  numerically  the  unsteady  small 
perturbation  transonic  equation.  In  this  approach,  a  set  of  locally  defined  trial 
functions  with  undetermined  parameters  is  assumed  as  the  appiroximate  solu¬ 
tion,  and  the  integral  expression  for  the  square  of  errors  committed  by  the 
approximate  solution  is  formulated.  Then  the  integral  of  square  errors  in 
the  time-space  domain  is  minimized  wtih  respect  to  the  undetermined  param¬ 
eters  to  yield  a  system  of  algebraic  equations.  In  actual  computations,  the 
minimization  process  is  performed  at  the  element  level  and  then  an  assembling 
process  is  invoked  to  obtain  the  system  of  algebraic  equations. 

Written  in  the  usual  manner  with  repeated  indices  implying  summation, 
the  approximate  solution  has  the  following  form 

<l>  -  N.  M.  <}k  (i  =  1  to  n,  k  =  1  to  3)  (11) 

1  K  1 
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where 


N.  =  N.  (x,y) 
1  i 


Mk=Mk(t) 


Herein,  N.  =  shape  functions  in  space,  Mk  =  shape  functions  in  time,  *,  =  the 
undetermined  parameter  at  node  i  and  time  level  k,  and  n  =  total  number 
of  unknown  parameters  at  one  time  level;  three  time  levels  are  involved 
in  the  time  marching  process.  As  discussed  in  the  first  volume  of  thus  re¬ 
port,  for  stability  reasons,  the  resulting  residual  is  modified  to  take  the 

following  form 

*  =  +  (U) 


where 


Br-^<2Mk,tNj,xtMk,ttNj) 


and  a  is  approximated  by 


a  =  (l-M*)-M*(l+r)N  *: 


p,x  p 


From  the  expression 
obtained  as 


of  R,  an  integral  expression  for  the  square  errors  is 


=  //Ri 


dAdt 


minimization  of  X  with  respect  to  the  undetermined  parameters  at  time 


Upon  minimization  oi  \ 
level  k  =  3,  one  obtains 


M-  R  dA  dt  =  0 

rv  A  * 


L  etting 


W  =  ~  =  (aN  +  N  )  -  B"  -  (1+  y)  N.  N  (19) 

1  i,xx  l,  yy  i  oo  l ,  x  I  ,  xx  I  '  ' 


OR 


3  ..2 


the  resulting  system  of  algebraic  equations  then  becomes 


f/w.  [<aNji 


«x  +  Nj,yy>-BJ-  ^dAdt  =  0 


(20) 


Or  in  the  form  of  a  recurrence  relationship,  as 


s  ..</>:  =  L. 

J  i 


where 


Sii  =  l/W  [ocN  +  N-  -  B31  dAdt 
J  J  1  1  e**  j.yy  J  J 


and 


(21) 


(22) 


Li  =  [“Nj.xx  +  Nj,yy  '  ^  I  ^  ^  dt,  k  =  2,3 


(23) 


Equation  (21)  will  be  used  to  solve  for 


3  12 

<t>.  in  terms  of  </>.  and  (/>.  . 
1  J  J 


As  stated  earlier,  the  system  matrix  is  obtained  by  combining  appropriate 
contributions  from  all  the  elements.  The  element  matrices,  in  turn,  are  eval¬ 
uated  effectively  by  numerical  integration  to  avoid  the  tedious  and  error  prone 
algebraic  manipulations.  The  Gaussian  quadrature  presently  used  has  seven 
points  in  the  spatial  direction,  and  two  points  in  the  time  direction.  This 


quadrature  formula  can  integrate  exactly  the  product  of  a  quintic  polynomial 
in  space  and  a  cubic  polynomial  in  time. 

1.  ELEMENT  DESCRIPTION 

The  basic  element  presently  used  is  the  product  of  an  element  in  space 
and  an  element  in  time.  For  the  finite  element  approximation  in  the  spatial 
directions,  the  nonconforming  cubic  triangular  elements  developed  by  Bazeley 
et  al.  (Ref.  5)  are  used.  Also  used  in  the  program  are  quadrilateral  elements 
constructed  from  these  triangular  elements.  The  element  in  time  is  a  quad¬ 
ratic  Lagrangian  element. 

The  basic  triangular  element  is  shown  in  Fig.  2,  which  at  each  vertex 
has  the  function  itself  and  its  two  first  derivatives  (velocity  components)  as 
undetermined  parameters.  This  type  of  element  was  adapted  mainly  because 
boundary  conditions  of  both  Dirichlet  and  Neumann  types  can  be  imposed  with 


Figure  2  -  Triangular  Element  with  Undetermined  Parameters 
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equal  convenience.  In  addition,  because  velocities  at  nodes  are  treated  as 
primary  unknowns,  secondary  unknowns  such  as  local  Mach  number,  and 
pressure  coefficients,  etc.,  can  be  computed  directly  without  resorting  to 
numerical  differentiation  and  thus  assuring  higher  accuracy.  Furthermore, 
the  use  of  higher  order  elements  can  usually  improve  computational  effi¬ 
ciency  as  evidenced  in  most  finite  element  analyses. 


In  the  element,  the  approximate  solution  is  assumed  as 


N.  0.  (i  =  1  to  9) 


(24) 


in  which  0  .'  s  are  the  nine  undetermined  parameters  of  <£  and  its  first  deriva¬ 
tives  at  the  nodal  points  as  shown  in  the  figure  and  N.'s  are  the  corresponding 
shape  functions,  which  are  expressed  in  terms  of  area  coordinates. 


Defined  in  the  following  are  these  shape  functions  and  their  first  and 
second  derivatives. 


Letting 


a.  =  x.y  -  x,  y . 

l  j7k  kyj 

b.  =  y.  -  y, 
i  J  yk 


c.  =  x,  -  x. 
1  k  J 


A  =  area  of  triangle  1-2-3  =  (b.ck  -  bkc  .)/2 

J  J 

a  =  0.5  (ck-c.) 

P  =  0.5  (b.-bk) 


H*  =  Vj^  +  V^i+Wj 

H*x  =  2^ibA +  Vkbi +  Wj> 
Hyy  ‘  2(hcjck  +  Vkci  +  5kcicj> 
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Nf  = 

h2<V 

k  -  Vj> 4 

■  PH 

Nf,x  = 

^|2bi^‘Vk- 

bk^  +  PHx] 

Ni.y  = 

Ta  l2' 

‘Wk- 

bkC  j)  +  2A;2  +  pHy 

Nl,xx  = 

(2A)2 

K'Vk  - 

-  Vj)  +  PHkx] 

pH  | 
P  YY, 

Ni,yy  = 

(2A)2 

K<Vk 

-  bk£j)  +4c.(2A);i  + 

The  undetermined  parameters  are  arranged  in  the  order  of  <^,  vp  U2’ 
v2’  ^3*  U3  and  V3‘ 

The  elements  constructed  from  the  basic  triangular  elements  are  shown 
in  Fig.  3  which  are  also  used  in  the  program.  The  quadrilaterals  are  used  to 
save  input  data  otherwise  necessary,  while  the  trapezoidals  can  be  used  to 
eliminate  bias  effects  and  to  remove  the  improper  upwind  influence  from  the 
downwind  station  in  the  calculation  of  steady  transonic  flow.  The  element 
matrix  for  a  quadrilateral  is  obtained  by  combining  approximately  the  matrices 
for  two  triangles,  while  the  matrix  for  a  trapezoidal  element  is  obtained  by 
combining  and  averaging  contributions  from  four  (two  left -running  and  two  right¬ 
running)  triangles.  In  any  event,  the  present  program  selects  automatically 
the  type  of  element  to  use,  according  to  the  flow  behavior  m  the  element. 

For  the  finite  element  approximation  in  time,  the  one  presently  used  is 
a  quadratic  element  using  Lagrangian  interpolation,  as  shown  in  Fig.  4.  In  the 
element,  the  approximate  solution  at  node  i  is  assumed  as 
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a.  Quadrilateral  Element 


b.  Trapezoidal  Element 


Elements  Constructed  from  Ba 


Figure  4  -  Shape  Functions  in  Time 


where  the  shape  function  M 


k 


is  defined  as 


m,  =  n 
k  Mi 


j  =  1  to  3,  k  =  1  to  3 


The  three  shape  functions  in  time  are  shown  schematically  in  Fig. 4. 


3.  THE  IMPOSITION  OF  BOUNDARY  CONDITIONS  AND  TIME  MARCHING 

As  stated  earlier,  the  imposition  of  boundary  conditions  for  the  present 
problem  can  be  carried  out  conveniently  because  the  elements  presently  used 
have  function  value  and  two  first  derivatives  as  primary  unknowns.  The  asso¬ 
ciated  boundary  conditions  thus  can  be  treated  as  the  essential  type,  i.e.,  having 
prescribed  values.  Standard  finite  element  methodology  is  therefore  followed 
by  assembling  first  an  unconstrained  problem  and  then  modifying  the  matrix 
equations  accordingly. 

On  the  airfoil,  the  boundary  condition  of  flow  tangency  is  imposed.  This 
is  accomplished  by  replacing  the  algebraic  equation  for  <j>  at  node  i  originally 

>  y 

generated  by 


(<t> 


y'i  = 


6  (h  + 


h,t>i,+ 


g  „  (1  +  $ 


x>i. 


(26) 


where  (</>  )j  designates  the  unknown 
condition  is  obtained  by  omitting  0 

9 

alo.  .-c  the  chordline. 


x 


at  node  "i".  A  linearized  boundary 
in  the  last  term,  which  is  to  be  imposed 


On  the  outer  boundary  of  the  finite  element  mesh,  for  computing  steady 
flow,  asymptotic  solution  in  the  form  of  Eq.(6)  is  imposed.  However,  as  it 
has  been  shown  that,  with  a  moderately  large  domain,  the  contribution  from 
the  thickness  and  the  area  integrals  are  of  higher  order  effects,  therefore 
only  the  first  term  representing  the  lifting  effect  is  incorporated  in  the 
present  program.  The  circulation  strength,  in  turn,  is  updated  systematically 


■ 


fli 
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according  to  the  potential  jump  calculated  at  the  trailing  edge.  For  unsteady 
flows,  the  circulation  strength  corresponding  to  mean  steady  flow  is  first 
computed  and  kept  unchanged  on  the  outer  boundary.  For  symmetric  airfoil 
without  angle  of  attack,  the  circulation  is  set  equal  to  zero. 


The  boundary  conditions  along  the  vortex  wake,  as  defined  by  Eqs.  (4a) 
and  (4b),  are  imposed  in  the  following  way.  To  impose  the  first  condition  for 
a  pair  of  nodes  along  the  wake,  the  two  finite  element  equations  generated 
originally  for  4>+  and  </>'  are  properly  combined  to  yield  one  equation,  while 

the  equation  for  <l>+  is  replaced  by  the  following  equation: 

t  X 


.(3)  (3) 

+  m3  t(</>  -  <t>  ) 


=  M. 


i.t<* 


(2)  ,(2) 

•  -  ) 


+  M 


L.t<*‘ 


(1)  ,(1) 

’  -  </>+  ) 


(27) 


where  M  is  the  shape  function  in  time  and  the  superscripts  denote  the  time 
steps,  while  solutions  for  the  first  and  the  second  time  steps  are  known. 

The  second  condition,  Eq.  (4b),  can  be  imposed  in  a  similar  fashion. 

This  condition,  within  the  limit  of  small  perturbation  assumption,  will  ensure 
that  flows  just  above  and  below  the  branch  cut  are  tangent  to  each  other.  There¬ 
fore,  flow  at  the  branch  cut  behaves  as  an  unsteady  vortex  sheet,  and  in  the 
steady  case  this  vortex  sheet  will  vanish  eventually. 

In  the  present  program,  no  special  efforts  have  been  devoted  to  treat 
the  leading  edge  singularity.  Any  rigorous  approach  should  consider  also  the 
invalidity  of  small  perturbation  theory  in  these  regions,  which  is  certainly 
beyond  the  scope  of  the  present  study.  Nevertheless,  these  regions  are  rela¬ 
tively  small  and  the  flow  is  usually  subsonic  in  nature,  which  implies  any  error 
committed  is  generally  localized  and  becomes  less  important  elsewhere.  A 
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possible  remedy  is  to  use  very  fine  elements  in  that  region  and  treat  the  lead¬ 
ing  edge  as  a  singular  joint.  In  practical  computations,  however,  an  element 
with  its  area  approaching  zero  might  create  a  new  numerical  problem  because 
all  shape  function  derivatives  involve  element  area  as  divisor.  Therefore, 
caution  must  be  taken  in  using  extremely  small  elements. 

With  the  equations  properly  assembled  and  boundary  conditions  imposed, 
the  system  of  algebraic  equations  is  finally  solved  by  marching  in  time  in  the 
form 


S..  <j> =  L. 
ij  J  i 


to  solve  for  the  solution  <f>  n  at  the  nth  time  step.  The  term  on  the  right-hand 
side  is  evaluated  according  to  Eq.  (23),  with  the  solutions  at  the  previous  two 
time  steps. 

When  a  steady  flow  is  to  he  computed,  Eq.(27)  is  to  be  solved  subject 
to  certain  prescribed  convergence  criteria.  The  one  presently  used  is  that 
the  relative  change  of  local  Mach  number  between  two  consecutive  time  steps 
should  be  less  than  a  prescribed  small  number  for  all  the  nodes  in  the  flow 
field,  that  is 


M(n>  -  M1"'1’ 


<  e 


Numerical  experimentations  indicate  that  a  solution  with  adequate  accuracy 
is  generally  obtainable  with  €  in  the  range  0.001  <  e  <  0.005. 


SECTION  IV 

PROGRAM  DESCRIPTIONS  AND  USAGE 

As  stated  earlier,  the  present  computer  program  is  capable  of  analyzing 
both  steady  and  unsteady  (oscillatory  and  transient)  transonic  flows  over  thin 
airfoils  by  solving  the  small  disturbance  but  nonlinear  transonic  potential  equa¬ 
tion.  The  numerical  algorithm  is  based  on  utilizing  the  finite  element  technique 
in  conjunction  with  the  least  squares  method  of  weighted  residuals.  To  solve 
a  particular  transonic  flow  problem,  an  appropriate  finite  element  mesh  must 
first  be  set  up,  followed  by  using  the  numerical  procedures  summarized  in 
Section  III.  Thus  the  present  programs  are  separated  into  two  parts  —  the 
first  part  which  generates  the  necessary  mesh  information  from  a  limited 
number  of  input  cards  and  the  second  part  which  carries  out  the  analysis  with 
element  mesh  generated  in  the  first  part.  By  doing  so,  the  generated  mesh 
can  be  fully  inspected  for  its  correctness  prior  to  the  analysis,  and  storage 
required  in  the  second  program  can  be  accordingly  determined  to  suit  each 
particular  problem.  These  two  programs  are  designated  UTRANL-I  and 
UTRANL-II  (Unsteady  Transonic  Flow  by  Lockheed,  parts  I  and  II,  respec¬ 
tively)  and  are  described  in  more  detail  in  the  following. 

1.  UTRANL-I 
Scope  and  Flow  Chart 

This  program  reads  in  a  limited  number  of  input  cards  and  generates 
mesh  information  such  as  element  nodes,  nodal  coordinates,  boundary  nodes 
and  airfoil  slope,  etc.,  to  be  used  in  the  second  program.  The  mesh  is  gene¬ 
rated  for  the  finite  flow  domain  indicated  in  Fig.  1,  with  the  outer  boundary 
set  to  be  a  circle  centered  at  the  midchord  of  the  airfoil.  The  present  pro¬ 
gram  generates  initially  a  mesh  based  on  the  geometry  of  a  6%  thick  circular 
arc  airfoil  and,  for  airfoils  of  other  geometric  shape,  the  y-coordinate  and 
surface  slope  of  nodes  on  the  airfoil  are  corrected  accordingly,  using  a  BLOCK 
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DATA  which  defines  the  geometric  shape  of  the  airfoil  under  consideration. 
Additionally,  the  present  program  has  an  option  for  plotting  the  generated 
mesh  for  quick  inspection.  A  schematic  flow  chart  of  the  program  is  shown 
in  Fig.  5. 

The  program  as  presently  dimensioned  requires  approximately  46g  K 
words  to  run  and  can  accommodate  the  following  maxima: 

•  400  elements 

•  400  nodal  points 

•  100  nodes  for  each  type  of  boundary  conditions,  and 

•  10  nodes  connected  to  any  node  in  the  mesh. 


Description  of  Variables 


TITLE 

IOPT 


NDEL 

X 

Y 

NIDS 

NID 

VAF 

NEM 

NPM 


Array  used  to  describe  the  problem 
under  consideration. 

Array  containing  the  option  keys  which 
are  activated  when  read  in  as  "  1."  Pre¬ 
sently  there  are  two  keys  in  the  program, 
IOPT  (9)  for  using  the  BLOCK  DATA 
to  generate  the  required  airfoil  geometry 
and  IOPT  (12)  for  mesh  plot  option. 

Array  containing  element  node  points. 

Array  containing  the  x-coordinate  of 
nodal  points. 

Array  containing  the  y-coordinate  of 
nodal  points. 

Array  containing  the  actual  number  of 
boundary  nodes  for  far  field,  on  the  line 
of  symmetry,  and  along  the  airfoil  surface. 

Array  containing  the  nodal  numbers  for 
each  type  of  boundary  nodes  mentioned 
above. 

Array  containing  the  airfoil  slope  for  nodes 
on  the  airfoil  surface. 

Maximum  number  of  elements  allowed. 
Maximum  number  of  nodal  points  allowed. 


Read 


Read 


START 


TITLE 

IOPT  (Control  Keys) 


MESHBW 


(Generate  mesh  infor¬ 
mation  from  a  limited 
number  of  input  cards, 
including  ID) 


NIDS,  NID 
(Boundary  nodes) 


Print  Output  Data 


Punch  Output  Data 
(or  store  on  tape) 


Back  to  Start  for  Next  Case' 
(Stop  if  EOF  read)  j 


* 

No 

Subroutine 

MESHBW 

Yes 

— ►<IOPT(9)  =  n>-^ 

Call 

RAIF 

Call  MPLOT 
(Plot  generated  mesh) 


Figure  5  -  Flow  Chart  of  UTRANL-I 


NELS 

NPS 

NBW 

ID 


XLE 

XTE 


Total  actual  number  of  elements. 

Total  actual  number  of  nodal  points. 

Full  bandwidth  of  the  resulting  system 
matrix. 

Control  key  for  calling  the  subroutine 
GOCR  to  compute  a  6%  thick  circular 
arc  airfoil  geometry  (y-coordinate  and 
surface  slope)  and  the  corresponding  far 
field  boundary  curve.  Values  on  the  upper 
surface  or  lower  surface  are  computed 
according  to  ID  =  1  or  -1.  GOCR  will  not 
be  called  if  ID  =  0. 

The  distance  between  the  leading  edge  and 
the  point  farthest  to  the  left  in  the  computa¬ 
tional  domain. 

The  distance  between  the  trailing  edge  and 
the  point  farthest  to  the  right  in  the  computa- 

rlnmain. 


Subroutines 

MESHBW:  This  subroutine  is  called  to  generate  mesh  information  in¬ 
cluding  the  following:  array  of  element  nodes,  total  number  of  elements.  x- 
and  y-coordinates  of  nodal  points,  and  the  maximum  difference  m  nodal  num¬ 
bers  between  two  connected  nodes.  In  generating  this  information,  a  small 
amount  of  input  cards  is  required,  which  will  be  described  later. 

MPBOT:  This  routine  plots  the  mesh  generated  together  with  the  original 
node  numbering,  using  the  CALCOMP  plotter.  Correctness  of  the  mesh  can 
thus  be  checked  quickly  prior  to  running  the  second  program. 

GOCR:  This  subroutine,  if  called,  generates  the  y-coordinate 
slop-  at  a  node  on  the  6%  thick  circular  arc  airfoil  by  specifying  the  corre¬ 
sponding  x-coordinate.  For  other  airfoils,  the  following  subroutines  and 
BLOCK  DATA  are  called  to  generate  the  correct  values  for  the  atrfotl  under 

consideration. 
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RAIF :  This  subroutine  is  called  when  IOPT(9)  =  1 ,  to  generate  the  y- 
coordinates  and  surface  slope  for  nodes  on  the  airfoil,  based  on  available 
tabulated  data  for  airfoil  thickness  distribution. 

GMPG:  This  subroutine  is  called  to  compute  the  surface  slope  for 
preselected  points  on  the  airfoil,  using  the  data  given  in  BLOCK  DATA. 

BLOCK  DATA:  This  data  block  defines  the  thickness  distribution  ac¬ 
cording  to  available  tabulated  data  for  the  specific  airfoil  to  be  studied.  This 
subroutine  is  to  be  used  for  any  airfoil,  provided  the  pertinent  data  statements 
are  supplied. 

Input  and  Output 

All  input  and  output  are  referred  to  the  automatic  numbered  scheme, 
with  input  in  the  form  of  punch  cards,  and  output  in  the  form  of  printout  and 
punch  cards.  If  subroutine  MPLOT  is  called,  the  mesh  is  also  plotted  for 
the  convenience  of  inspection. 

Input 


Input  cards  to  this  program  should  be  prepared  and  provided  in  the 
order  described  below. 


A.  TITLE  CARD  (12A6) 

Col.  1-72  Description  of  the  problem  under  study 

B.  OPTIONS  CARD  (4012) 

Col.  18  Punch  "  1"  for  airfoil  other  than  the  6%  thick 
circular  arc. 

Col.  24  Punch  "  1"  if  mesh  plot  is  desired 

C.  ELEMENT  CARDS  (1615) 

These  are  cards  supplied  to  generate  element  nodal  numbers 
in  groups.  One  card  is  needed  for  each  group  of  elements 
whose  nodes  are  related  in  a  regular  pattern.  The  number 
of  elements  in  a  group  can  vary  from  one  to  any  positive 
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integer  number,  depending  on  how  the  nodes  are  numbered 
originally.  There  can  be  as  many  cards  as  required  to 
generate  the  element  nodes,  and  a  blank  card  is  added  at 
the  end  to  terminate  the  process. 


I- 5  Number  of  nodes  per  element.  Punch  3  for 

triangular  elements,  and  punch  4  for  quadri¬ 
lateral  and  trapezoidal  elements. 

6-10  Number  of  elements  in  one  direction  (called 

first  direction). 

II- 15  Number  of  elements  in  the  other  direction 

(called  second  direction). 

16-20  Increment  of  nodal  numbers  in  the  first 
direction. 

21-25  Increment  of  nodal  numbers  in  the  second 
direction. 

26-30  First  nodal  number  in  the  element 

31-35  Second  nodal  number  in  the  element 

36-40  Third  nodal  number  in  the  element 

41-45  Fourth  nodal  number  in  the  element 
("0"  for  triangles) 


As  stated  before,  element  nodes  are  presently  ordered  in  the  counterclock¬ 
wise  direction,  and  the  upper  left  corner  node  should  be  taken  as  the  first 
nodal  point  in  an  element  located  in  the  anticipated  supersonic  region.  This 
is  required  to  enact  the  assembling  process  correctly  in  the  supersonic  pocket 
for  steady  flow. 

For  example,  to  generate  element  nodes  for  the  mesh  shown  on  page  25, 
three  input  cards  are  required:  one  card  for  the  first  12  elements,  second  card 
for  the  next  two  elements,  and  the  third  card  for  the  triangular  element.  The 
three  cards  are; 

44341  2  1  5  6 

4  1  2  0  1  18  17  21  22 

3  1  1  0  0  20  19  23  0 

(Ended  by  a  blank  card  when  all  elements  are  covered) 
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D.  CARD  FOR  TOTAL  NO.  OF  NODES  (15) 

Col.  1-5  Total  number  of  nodes  for  the  entire  flow  region 

E.  CARD  FOR  LEADING  AND  TRAILING  EDGE  (2F10.0) 

Col.  1-10  The  distance  between  the  leading  edge  and  the 
point  farthest  to  the  left  in  the  computational 
domain. 

11-20  The  distance  between  the  trailing  edge  and  the 
point  farthest  to  the  right  in  the  computational 
domain. 

F.  NODE  COORDINATE  CARDS  (315,  5F10.0,  15) 

The  x-  and  y-coordinates  of  nodes  are  also  generated  in 
connected  groups.  A  connected  group  of  nodes  are  those 
falling  on  a  straight  line  with  a  constant  ratio  of  nodal  dis¬ 
tances  and  constant  increment  in  nodal  numbers  as  well. 

The  intermediate  nodal  points  can  thus  be  generated  by 
linear  interpolation.  Again,  the  process  is  terminated 
by  encountering  a  blank  card. 

Col.  Description 

I- 5  Node  number  of  the  beginning  nodal  point. 

6-10  Total  number  of  nodes  on  the  line. 

II- 15  Increment  of  nodal  numbers  between  two 

consecutive  nodes. 


16-25  Ratio  of  distances  between  three  consecutive  nodes. 
26-35  x-coordinate  of  the  beginning  node. 

36-45  y-coordinate  of  the  beginning  node. 

46-55  x-coordinate  of  the  end  node. 

56-65  y-coordinate  of  the  end  node. 

66-70  0  if  a  node  does  not  fall  on  the  airfoil  surface. 

1  for  calling  GOCR  to  compute  y-coordinate  and 
surface  slope  for  the  upper  surface;  -1  for  those 
on  the  lower  surface. 

For  example,  to  generate  the  coordinates  for  the  nodes  depicted  below  with 
a  constant  distance  ratio  of  1.2,  the  input  card  should  read  as 


5  4  1.2  y  z  X 18  y18 


G.  CARD  FOR  NO.  OF  NODES  ON  BOUNDARIES  (1615) 
Col.  Description 

I- 5  Number  of  nodes  on  the  outer  boundary 

(called  NFARF) 

6-10  Number  of  nodes  on  line  of  symmetry 
(called  NWAKE) 

II- 15  Number  of  nodes  on  the  airfoil  surface 

(called  NBODY) 


H.  CARDS  FOR  BOUNDARY  NODES  (1615) 

(a)  Nodes  on  the  outer  boundary  (NFARF  entries) 

(b)  Nodes  on  line  of  symmetry  (NWAKE  entries) 


Output 

The  output  from  this  program  is  in  the  form  of  printout,  punch  cards, 
and  also  a  plot  if  the  mesh  plot  option  is  invoked. 

The  following  items  are  printed  in  the  order  as  described: 

•  Title  of  the  problem  under  study 

•  Control  keys  specified  in  the  options  card 

•  Total  number  of  elements,  number  of  nodal 
points,  and  the  full  bandwidth 

•  Element  numbers  and  element  node  points 

•  Nodal  numbers  and  their  corresponding  x- 
and  y -coordinates 

•  Nodes  on  each  type  of  boundaries,  and 

•  Slope  for  nodes  on  the  airfoil  surface. 

The  above  items  ,  except  the  first  two,  are  then  punched  and  saved  as  input  to 
the  second  program  (UTRANL-II). 

2.  UTRANL-II 
Scope  and  Flow  Chart 

This  program  carries  out  the  analysis  following  the  numerical  pro¬ 
cedures  described  in  Section  III,  with  any  mesh  information  generated  and 
supplied  through  the  first  program.  In  summary  this  program  solves  the 
small  perturbation  transonic  potential  equation,  by  marching  in  time,  via 
the  least  squares  finite  element  technique.  The  embedded  shock  waves  are 
smeared  and  the  resulting  system  of  algebraic  equations,  relating  the  solu¬ 
tion  of  current  time  step  to  the  solutions  of  previous  two  time  steps,  is  solved 
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by  direct  elimination.  The  present  program  has  four  options,  one  for  inputting 
nonzero  initial  guess  of  the  solutions,  one  for  computing  in  the  same  run  solu¬ 
tions  for  higher  Mach  number(s)  using  results  computed  for  lower  Mach  number 
case(s),  the  third  one  for  selecting  the  form  of  boundary  conditions  on  the  air¬ 
foil,  and  the  last  one  for  updating  the  circulation  strength.  A  schematic  flow 
chart  of  the  program  is  shown  in  Fig.  6. 

The  program  as  presently  dimensioned  requires  approximately  177gK 
words  for  the  program  itself  and  can  accommodate  the  following  maxima: 

•  173  elements 

•  202  nodal  points 

•  50  nodes  for  each  type  of  boundary  conditions,  and 

•  84  in  full  bandwidth  of  the  resulting  matrix 
equations . 

Description  of  Variables 


TITLE  Array  used  to  describe  the  problem  under 
consideration 

IOPT  Array  containing  the  option  keys  which  are 

activated  when  read  in  as  "  1."  Presently 
there  are  four  options:  IOPT(l)  for  con¬ 
tinuation  for  high  Mach  number  cases  while 
using  existing  results  computed  for  earlier 
case,  IOPT(2)  for  inputting  non-zero  initial 
guess,  IOPT(3)  for  selecting  the  linearized 
boundary  condition  applied  along  the  airfoil 
chordline,  and  IOPT(4)  for  updating  the 
circulation  strength. 

NOD  Array  containing  element  node  points 

S  The  resulting  system  matrix  with  main 

diagonal  terms  stored  in  the  column 
numbered  NHBW.  After  solving,  the 
column  numbered  NBW  contains  the 
solution  for  the  current  time  step. 

SLl  Array  containing  the  solution  for  the  first 

time  step 

SLP  Array  containing  the  solution  for  the  second 

time  step. 
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Figure  6  -  Flow  Chart  of  UTRANL-I1  (Completed) 
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Array  containing  the  x-coordinate  of  nodal 
points 

Array  containing  the  y-coordinate  of  nodal 
points . 

Array  containing  values  of  local  Mach  number 
at  nodal  points  for  the  current  time  step. 

2  2 

Array  containing  values  of  +  M^  (1  +  7)u 
evaluated  at  nodal  points. 

Array  containing  the  actual  number  of  boundary 
nodes  for  far  field,  on  the  line  of  symmetry, 
and  along  the  airfoil,  respectively. 

Array  containing  the  nodal  numbers  for  each 
type  of  boundary  nodes  mentioned  above. 

Array  containing  the  airfoil  slope  for  nodes 
on  the  airfoil  surface. 

Array  to  store  values  for  use  in  subroutine 
FPLOT. 

Another  array  to  be  used  in  the  subroutine 
FPLOT. 

Maximum  number  of  elements  allowed. 

Maximum  number  of  nodal  points  allowed. 

Maximum  full  bandwidth  allowed  in  the 
resulting  system  matrix. 

Maximum  number  of  equations  allowed. 

Total  actual  number  of  elements. 

Total  actual  number  of  nodal  points. 

Full  bandwidth  of  the  resulting  system 
matrix. 

Half  bandwidth  of  the  resulting  system 
matrix. 

Number  of  equations  in  the  resulting 
system  of  equations. 

Small  number  used  to  check  convergence 
based  on  the  relative  change  of  local  Mach 
number  between  two  consecutive  iterations 
in  steady  flow  computations. 

Freestream  Mach  number. 
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RKSC 


Reduced  frequency  based  on  semi-chord, 
and  then  converted  into  that  based  on  chord 
length  in  computations.  (RKSC  =  0  for  steady 
flow). 

ANDG  Amplitude  in  degrees  of  the  oscillating  flap. 

(ANDG  =  0  for  steady  flow). 

XHP  x-coordinate  of  the  hinge  point. 

XLE  x-coordinate  of  the  trailing  edge. 

DTM  Desired  number  of  time  steps  in  each  period 

of  time,  which  is  then  converted  into  time 
increment  At  in  computations.  (For  steady 
flow  the  input  value  of  DTM  is  used  as  At). 

TFIN  Total  number  of  time  steps  to  be  computed. 

TO  Initial  time  level,  usually  set  equal  to  zero. 

ALPHA  Angle  of  attack  of  the  airfoil  at  mean  steady 
position. 

CIRCO  Circulation  strength  for  airfoil  at  its  mean 
position. 

SQMAC  Square  of  the  freestream  Mach  number. 

IRES  Number  of  time -marching  steps. 

POT  Perturbed  velocity  potential. 

UPT  Perturbed  velocity  component  in  the  x- 

direction. 

V  Perturbed  velocity  component  in  the  y- 

direction. 

U  Total  velocity  component  in  the  x-direction. 

PRA  The  ratio  of  local  static  pressure  to  stagnation 

pressure. 

CP  Pressure  coefficient. 

DELM  Change  of  local  Mach  number  between  two 

consecutive  time  steps. 

CPU  Normalized  unsteady  pressure  coefficient. 

Subroutines 

NEWK:  This  subroutine  is  called  in  the  main  program  to  generate  system 
matrix  for  the  entire  flow  field  by  assembling  appropriately  the  element  matrices. 
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The  elements  are,  in  general,  a  combination  of  triangles,  quadrilaterals  and 
trapezoids.  For  a  4-node  element,  this  subroutine  also  determines  either  to 
treat  it  as  a  quadrilateral  or  as  a  trapezoid,  depending  on  the  local  behavior 
of  the  governing  equation  being  elliptic  or  otherwise.  Subroutine  EMQT  is 
called  to  gene  rate  the  element  matrices. 

OUTP:  This  subroutine  is  called  in  the  main  program  to  compute  the 

secondary  parameters  (Mach  number,  pressure  coefficient,  etc.),  and  to  print 

all  the  computed  results  at  all  nodes.  These  include  the  perturbed  velocity 

potential,  total  velocity  components  in  the  x-  and  y-directions ,  the  value  of 
2 

M  (1  +  2.4u),  local  Mach  number,  the  ratio  of  local  static  pressure  to  stag- 
oo 

nation  pressure  coefficient,  and  finally  the  change  of  local  Mach  number 
between  two  consecutive  time  steps.  Also  computed  are  the  differences  of 
pressure  coefficients  between  the  upper  and  lower  airfoil  surfaces. 

BCDS:  This  subroutine  is  called  in  the  main  program  to  impose  all  the 
relevent  boundary  conditions.  These  conditions  include  the  far  field  condition, 
the  airfoil  flow  tangency  condition,  and  the  condition  of  equal  pressure  along 
the  wake. 

EMQT ;  This  subroutine  generates  element  matrix  for  elements  in  the 
form  of  triangles,  quadrilaterals  and  trapezoids.  Element  matrix  for  the 
latter  two  types  is  obtained  by  combining  appropriately  (and  averaging  for  the 
last  type)  contributions  from  triangular  elements  which  are  in  turn  evaluated 
in  subroutine  EMTC. 

EMTC:  This  subroutine  evaluates  the  element  matrix  for  a  space-time 
element  by  numerical  integration,  presently  with  Gaussian  quadrature.  The 
one  used  in  the  program  has  7  points  in  space  and  2  points  in  time,  which  inte¬ 
grates  exactly  a  quintic  function  in  space  and  a  cubic  function  in  time. 
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Shown  below  are  the  locations  of  the  Gaussian  points  in  space  with  their 
corresponding  weights. 


Points  Area  Coordinates 


1/3  1/3  1/3 

Pi  Pi 


m 


al  =  0.05961587 


=  0.47014206 


a2  =  0.79742699 


02  =  0.10128651 


Weights 


0.225 


0.13239415 


0.12593918 


Shown  below  are  the  Gaussian  points  in  time  with  their  corresponding 
weights. 


-1 

p  0  q 

1 

Points 

Normalized  Coordinate 

Weights 

P 

-0.55555556 

0.5 

q 

+0.55555556 

0.5 

TDRU;  This  is  a  subroutine  called  by  EMTC  to  evaluate  the  first  and 
second  derivatives  of  the  shape  functions  in  time  at  a  Gaussian  point. 

DERV:  This  is  a  subroutine  called  by  EMTC  to  evaluate  the  first  and 
second  derivatives  of  the  shape  functions  in  space  at  a  Gaussian  point,  based 
on  the  expressions  listed  in  Section  III  (2). 

BNDEQ:  This  is  an  equation  solver  for  banded  nonsymmetric  system 
of  algebraic  equations,  utilizing  the  direct  Gaussian  elimination  scheme.  In 
this  solver  the  system  matrix  is  arranged  in  a  twisted  form  such  that  main 
diagonal  terms  are  stored  in  the  column  numbered  NHBW  and  the  right-hand 
side  vector  is  in  the  column  NBW.  After  solving,  the  column  numbered  NBW 
contains  the  solution  vector. 

FP LOT :  FPLOT  collects  and  stores  data  as  it  becomes  available,  and 
upon  signal,  produces  a  printer  plot  in  practically  any  orientation  and  size. 

It  should  be  regarded  as  a  general  purpose  output  routine  for  displaying  out¬ 
put  data  in  graphical  form. 
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Description  of  Parameters 


Ml 

IPNT 

AR 

LR 

ISTOP 

NC 

NCMAX 

VI 

V2 


Size  of  the  main  storage  array,  and  it 
should  never  be  larger  than  200,  which 
corresponds  to  100  points  to  be  plotted. 

Counter  initialized  —  usually  IPNT  =  0  - 
in  the  calling  program.  It  is  incremented 
by  2  each  time  a  new  data  point  is  entered 
in  AR. 

Main  storage  array.  It  should  be  in  a 
dimension  statement  in  the  calling  pro¬ 
gram.  For  example,  DIMENSION  AR(800). 

An  array  of  bytes  used  to  hold  the  curve 
number.  It  should  be  dimensioned  for  Ml/2. 
The  type  declaration  LOGICAL*l  LR(400) 
should  be  in  the  calling  program. 

Flag  used  to  signal  that  all  data  has  been 
entered.  ISTOP  =  0  causes  data  to  be 
stored.  If  ISTOP  =  -1  and  NC  =  NCMAX 
the  program  immediately  branches  to  the 
plotting  section.  If  ISTOP  =1  and  NC'  = 
NCMAX  a  data  point  is  stored  and  then  the 
plotting  section  is  entered. 

Curve  number  for  which  data  is  being  en¬ 
tered.  It  must  be  a  positive  integer  less 
than  21. 

The  number  of  curves  to  appear  on  the  graph. 
This  is  the  largest  value  NC  will  have. 

The  horizontal  coordinate  of  the  data  point  to 
be  plotted. 

The  vertical  coordinate  of  the  data  point  to 
be  plotted. 


Input  and  Output 


Input 

The  bulk  of  input  to  this  program  is  provided  through  UTRANL-I  in  the 
form  of  punch  cards.  To  run  several  cases  of  various  freestream  Mach  num¬ 
ber  and  oscillating  frequency  but  using  the  same  mesh  and  boundary  conditions, 
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four  additional  cards  should  be  furnished  for  each  case  to  be  computed.  These 
cards  should  be  provided  in  the  following  order 

A.  TITLE  CARD  (12A6) 

Col.  1-72  Description  of  the  problem  under  study 

B.  OPTIONS  CARD  (4012) 

Col.  2  Punch  "1"  if  it  is  a  continuation 
in  computation  from  one  case 
to  another. 

Col.  4  Punch  "1"  if  cards  for  nonzero  initial 
guess  is  furnished  in  the  input  data. 

If  so,  insert  them  following  those 
generated  from  UTRANL-I. 

Col.  6  Punch  "1"  if  linearized  boundary  condi¬ 
tion  along  the  chordline  is  desired. 

Col.  8  Punch  "  1"  if  circulation  is  to  be  updated. 

NOTE:  As  seen  in  the  flow  ohart,  when  IOPT(l)  =  1,  the  next  two  options 
become  inoperative. 

C.  PROBLEM  PARAMETERS  CARDS  (8F10.0) 

First  Card 

Col.  1-10  Freestream  Mach  number. 

Col.  11-20  Reduced  frequency  based  on  semi-chord 
(0  for  steady  flow). 

Col.  21-30  Amplitude  in  degrees  for  the  oscillating 
flap.  (0  for  steady  flow). 

Col.  31-40  x-coordinate  of  the  hinge  point. 

Col.  41-50  x-coordinate  of  the  trailing  edge. 

Col.  51-60  Number  of  time  steps  desired  in  one  period 
of  time  or  time  step  size  for  steady  flow. 

Col.  61-70  Total  number  of  time  steps  to  be  computed. 

Col.  71-80  Initial  time  level  counted  in  time  steps. 

(Usually  set  equal  to  zero;  with  nonzero  value 
when  initial  solutions  are  provided) 


Second  Card 

Col.  1-10  Small  number  used  to  check  conver¬ 
gence  based  on  the  relative  change  of 
local  Mach  number  between  two  con¬ 
secutive  time  steps,  if  flow  approaches 
a  steady  flow,  usually  taken  in  the 
range  from  0.005  to  0.001. 

Col.  11-20  Mean  angle  of  attack  of  the  airfoil 

Col.  21-30  Circulation  for  flow  with  the  airfoil 
placed  at  its  mean  position. 


For  the  first  case  of  the  run,  punch  cards  from  UTRANL-I  are  placed  after 
the  above  cards.  For  subsequent  cases  using  the  same  mesh,  four  additional 
cards  as  described  above  are  required  for  each  case. 


Output 

The  output  from  this  program  include  printout,  punch  cards  and  solutions 
for  the  last  two  time  steps  (or  convergent  solutions  for  steady  flow). 


The  printouts  are  in  the  following  order: 


•  Title  of  the  problem  under  study 

•  Convergent  limit  (for  steady  flow) 

•  Control  keys  specified  in  the  options  card 

•  Total  number  of  element,  number  of  nodal  points, 
and  the  full  bandwidth  of  the  resulting  matrix 
equations 

•  Element  numbers  and  element  node  points 

•  Nodal  numbers  and  their  corresponding  x-  and  y- 
coordinates 

•  Boundary  nodes  for  each  type  of  boundaries 

•  Slope  for  nodes  on  the  airfoil  surface 

•  For  each  time  step,  the  computed  results  at  all 
nodal  points  are  printed.  These  include  the  per¬ 
turbed  velocity  potential,  total  velocity  components 

in  the  x-  and  y-directions,  the  value  of  M'j  (l  +  2.4u). 
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local  Mach  number,  the  ratio  of  local  static  pressure 
to  stagnation  pressure  coefficient,  and  finally  the 
change  of  local  Mach  number  between  two  consecutive 
time  steps. 

•  For  each  time  step,  the  unsteady  pressures  on  the  air¬ 
foil,  for  the  upper  and  lower  surfaces,  as  well  as  their 
difference,  are  printed.  At  the  end  of  each  period,  the 
mean  pressures  on  the  airfoil,  for  the  upper  and  lower 
surfaces,  obtained  by  time  average  within  a  period,  and 
their  difference,  are  also  printed. 

•  For  steady  case,  the  value  of  Cp  along  the  airfoil  is  dis¬ 
played  graphically.  For  unsteady  case,  the  difference 
of  unsteady  pressures  across  the  airfoil  is  displaced 
graphically. 


SECTION  V 
SAMPLE  CASES 


Two  sample  cases  are  presented  herein  to  demonstrate  the  usage  of 
the  computer  programs  described  in  the  previous  section  regarding  input 
and  typical  output  from  the  programs.  The  two  problems  are  flow  over  a 
NACA  64  A006  airfoil  with  a  quarter -chord  oscillating  flap  and  flow  over  a 
lifting  NACA  64  A410  airfoil. 

In  computing  the  steady  or  unsteady  transonic  flow  field  for  a  given 
airfoil  shape  with  the  present  programs,  the  following  procedures  are  gener¬ 
ally  followed: 


•  A  desired  mesh  is  sketched  with  each  node  assigned  a  number. 
In  order  to  save  storage  and  computation  time,  nodal  points 
should  be  numbered  along  the  direction  with  less  number  of 
nodes,  such  as  shown  in  Fig.  7. 

•  Based  on  the  mesh  sketched,  an  appropriate  set  of  input  cards 
is  prepared  and  supplied  to  UTRANL-I  to  generate  the  required 
mesh  information  in  the  form  of  punch  cards. 

•  The  above  punch  cards,  with  four  additional  cards  for  each 
case  and  cards  for  non-zero  initial  guess  if  available,  are 
used  as  input  to  UTRANL-II  to  compute  the  flow  field. 

•  During  the  execution  of  UTRANL-II,  results  for  each  time 
step  are  printed,  and  solution  for  the  last  two  time  steps 
is  saved  in  the  form  of  punch  cards  for  possible  later  use. 


The  following  paragraphs  describe  in  more  detail  the  input  and  output 
for  the  two  sample  problems. 

1.  FLOW  OVER  A  NACA  64  A006  AIRFOIL  WITH  AN  OSCILLATING  FLAP 

This  problem  was  analyzed  using  the  mesh  shown  in  Fig.  7,  which  con¬ 
sists  of  173  elements  together  with  202  nodes.  Summarized  below  are  the 
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input  and  output  to  theUTRANL-I  and  UTRANL-U  programs  for  the  present 
problem. 


UTRANL-I 


Input 

The  necessary  set  of  input  cards  is  listed  on  the  next  two  pages.  Note 
that  both  options  for  plotting  mesh  and  using  BLOCK  DATA  are  to  be  enacted, 
and  the  final  mesh  will  have  nodes  on  the  actual  airfoil  surface  rather  than 


on  the  chordline. 
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Input  Cards  to  UTRANL-I  for  NACA  64A006  Airfoil 

UNbTtADy  TRANSONIC  F LOW — MESH  t_3  MUOIFItU — 173  hUcMtNTS  —  202  NODtS 
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14 
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4 

1 
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5 
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3 
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2 
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4 
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14 

7 
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33 
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1 
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37 
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As  stated  before,  output  of  this  program  includes  printouts,  punch  cards 
and  a  plot  for  the  generated  mesh. 

Printouts  and  punch  cards  for  the  present  mesh  are  listed  on  the  next 
11  pages. 
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• 23733E+0 1 

173 

•12600E+0 1 

- • 2 1600  E  +  3 1 

174 

•93U90E+0 0 

- . 11496E+0 1 

175 

•  743  66E  +  Q  0 

- . 55528E+00 

176 

•  63117E  +  0  0 

-.2C56EE«-00 

177 

•56500E+0 0 

-.1C00QE-04 

178 

•  5650  0E  +  0  0 

.10000E-04 

179 

.63U7E  +  00 

•  2C566E  +  0  0 

180 

•  74366E+0  0 

. 55528E+Q0 

181 

•93490E+00 

•  1 1496E  +  0 1 

182 

•  12600E  +  0  1 

•  216CC  E+01 

183 

.17  300E+0  1 

- • 1 8C0C  E  +  0 1 

184 

• 11949E+0 1 

-.  86941E+00 

185 

.88Q15E+0  0 

- • 322C1E+00 

186 

. 6950  0E  +  0  0 

-.1CC00E-04 

187 

.6950 OE+O 0 

. 10C00E-04 

188 

.88015E+00 

. 32201E+00 

189 

• 11949E+0 1 

• 86941E+00 

190 

.17300E+01 

.1800CE+01 

191 

.  21400E*0  1 

-.130CCE+01 

192 

•  1 378  2E  +  0 1 

-.46429E+0Q 

193 

.9550QE+00 

- • 1 G0GQE-Q4 

194 

. 95500E  +  Q  0 

.15000E-04 

195 

•13782E+01 

.46429E+00 

196 

.  2140  0E  +  0 1 

.13C00E+01 

197 

.2380QE+0  ! 

-.78CCCE+00 

198 

•  16  0  OCE  +  O  1 

-.1Q00CE-34 

199 

.  160  0  0E  +  0  1 

. 10C00E-04 

200 

•  23800E  +  0  1 

•  7  8C00E  +  00 

201 

•  25000E+0 1 

-.1C000E-04 

202 

,  25  0  0  0E  +  0  1 

• 1COO0E-O4 
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173 

202 

84 

42 

10 

22 

4 

1 

2 

3 

3 

2 

7 

8 

24 

23 

34 

35 

2 

1 

5 

6 

7 

6 

12 

13 

14 

7 

13 

0 

14 

13 

21 

22 

23 

14 

22 

0 

22 

21 

31 

32 

23 

22 

32 

33 

2b 

24 

36 

37 

26 

2b 

37 

38 

lb 

24 

2b 

0 

16 

15 

25 

26 

8 

15 

16 

0 

9 

8 

16 

17 

4 

3 

9 

10 

30 

29 

41 

42 

33 

32 

44 

45 

34 

33 

45 

46 

44 

43 

55 

56 

4b 

44 

56 

57 

85 

54 

66 

67 

56 

55 

67 

68 

66 

65 

77 

78 

67 

66 

78 

79 

70 

69 

81 

82 

7b 

77 

89 

90 

81 

80 

92 

93 

82 

81 

93 

94 

92 

91 

103 

104 

93 

92 

104 

105 

103 

102 

1  14 

115 

104 

103 

1  15 

1  16 

114 

113 

12b 

126 

115 

1  14 

126 

127 

118 

117 

12* 

130 

l<£b 

12b 

137 

136 

129 

128 

140 

141 

130 

129 

141 

142 

140 

139 

151 

152 

141 

140 

152 

153 

lb  1 

150 

162 

163 

152 

151 

163 

164 

36 

35 

47 

48 

37 

36 

48 

49 

40 

39 

51 

52 

48 

47 

59 

60 

51 

50 

62 

63 

52 

51 

63 

64 

62 

61 

73 

74 

63 

62 

74 

7b 

73 

72 

84 

85 

74 

73 

85 

86 

84 

83 

95 

96 

85 

8*» 

96 

97 

88 

87 

99 

100 

96 

95 

107 

106 

99 

98 

110 

111 

100 

99 

1  1  1 

1  12 

1  10 

109 

121 

122 

111 

1  10 

122 

123 

121 

120 

132 

133 

122 

121 

133 

134 

132 

131 

143 

144 

133 

132 

144 

14b 

136 

135 

147 

148 

144 

143 

155 

156 

147 

146 

158 

159 

148 

147 

159 

160 

158 

157 

169 

170 

159 

158 

170 

171 

163 

162 

174 

1 7b 

164 

163 

17b 

176 

174 

173 

183 

184 

17b 

174 

184 

185 

184 

183 

191 

192 

lab 

184 

192 

193 

193 

192 

198 

0 

198 

197 

201 

0 

170 

169 

179 

180 

171 

170 

180 

181 

180 

179 

187 

188 

161 

180 

188 

169 

189 

188 

194 

1  95 

190 

189 

195 

196 

200 

199 

202 

0 

2»440E+00-5«270E-Ol-i«*20E+00-2.550E-0 l* 
2«035E+00-l •460E+00-i«263E+00-5»7682-0 l- 
1.283E+00  5.788E-01-*»035E+00  1.460E+00' 
a* 164E-0 1-3«455E-Oi-6»450E-Ol-4.500E-02- 

1 • 1 2SE+00  8»a64E— 0 i~i »6B0E+00  1 •660E+00' 
6»B80E-01“5»^33E-Oi-5«917E-Oi-2»323E-0 
5«9 17E— o 1  2  »323E-0 1—8 •68oE— 0 i  5.933E-01- 

7.7o5E-ol-2.37dE+00-6»4b2E-Ol-l»328t+00- 
5.049E-0  1-1*  o24t-0i-4.y00E-0l-t>*620E-03' 
5»  J02E-0 1  3»462E-0 i— 5  «732E— 0 1  7»096c-0l- 
7.o69E-ol-2»3y7E+00-i>»942e-ol-l»342t+oO- 
4.O37E-01-1 .4o4t-0i-4.500E-Ol-l »399E-02 
4.B71E-Q1  3.5bJE-0i-S«2b7E-0l  7.2Q5E-01 


a 

7 

14 

15 

15 

14 

23 

24 

7 

2 

6 

0 

6 

5 

1  1 

12 

12 

1  1 

19 

20 

13 

12 

20 

21 

20 

19 

29 

30 

21 

20 

30 

31 

34 

23 

33 

0 

24 

35 

36 

0 

27 

26 

36 

39 

28 

27 

39 

40 

17 

16 

26 

27 

lo 

17 

27 

26 

10 

9 

17 

18 

3 

6 

9 

0 

31 

30 

42 

43 

32 

31 

43 

44 

42 

41 

53 

54 

43 

42 

54 

55 

46 

45 

57 

58 

54 

53 

65 

66 

57 

56 

66 

69 

58 

57 

69 

70 

68 

67 

79 

80 

69 

68 

80 

81 

79 

78 

90 

91 

80 

79 

91 

92 

90 

89 

101 

102 

91 

90 

102 

103 

94 

93 

105 

106 

102 

101 

113 

1  14 

10b 

104 

1  16 

117 

106 

105 

1  17 

1  18 

1  16 

1  15 

127 

126 

1  17 

1  16 

128 

129 

127 

126 

1  86 

139 

128 

127 

139 

140 

136 

1  j7 

149 

150 

139 

136 

150 

151 

142 

141 

153 

154 

150 

149 

161 

162 

153 

152 

164 

165 

154 

153 

165 

166 

36 

37 

49 

50 

39 

38 

50 

51 

49 

48 

60 

61 

50 

49 

61 

62 

60 

59 

71 

72 

61 

60 

72 

73 

64 

63 

75 

76 

72 

71 

83 

84 

75 

74 

86 

87 

76 

75 

87 

88 

86 

85 

97 

98 

87 

86 

98 

99 

97 

96 

106 

109 

9b 

97 

109 

1  10 

108 

107 

1  19 

120 

109 

108 

120 

121 

112 

1  1  1 

123 

124 

120 

119 

131 

132 

123 

122 

134 

135 

124 

123 

135 

136 

134 

133 

145 

146 

135 

134 

146 

147 

14b 

144 

156 

157 

14b 

145 

157 

158 

156 

155 

167 

168 

157 

156 

168 

169 

160 

159 

171 

172 

162 

161 

173 

174 

16b 

164 

176 

177 

166 

165 

177 

0 

176 

175 

18b 

186 

177 

176 

186 

0 

186 

16b 

193 

0 

192 

191 

197 

198 

168 

167 

178 

0 

169 

166 

178 

179 

172 

171 

161 

182 

179 

176 

187 

0 

162 

181 

189 

190 

166 

187 

194 

0 

195 

194 

199 

0 

196 

195 

199 

200 

1  * 420E+00  2»850E-01-2»440E+00  5»270E-01 
8.870E-01-1 • 150E— 01 -8»670E-0 1  1 • 150E-0 1 
■  1 •680E+00”! »860E+00~l • l25E-*-00-6»864E-0 1 
■6»450E-01  4»500E-02-8» 164E-0 1  3»455E-01 
■1 • 130E+00-2»2SOE+00-8»517E-01-1 »207E+00 
■5»3S0E-0l“2»000E-02-5»350E-0 1  2»000E-02 
■8.517E-01  1  »207E +00-1  •  130E>00  2«250E+00 
-5»732E-0 l-7»098L-0 l-5«302E-0 l-3»462E-0 1 
■4  «900E-0 1  O»620E-03-5«049E-0 1  1»324E-01 
■&.462E-01  1  »328E4-00-7»705E-0l  2.378E4-00 
-b»267E-0 l-7»20bE-0l-*«87 IE-01 -3.553E-01 
-4.800E-01  l  .399L-02-4»t»37E-0t  1.40*E-01 
-S.942E-01  l  t342E4-00-7»0o9E-0 1  2»397E+00 
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5.yEbE-ol-2.42y8+00-4.9o2E-0l-l«38EE+00- 
3.obbE-o 1-1 .487E-0A-3.75OE-0  »-2.  101E-02- 

4  •  Ob  1 E— o  1  j.O=bE-0i-^»->‘>5E-0  1  7.3462-01- 
3.y658-ol-2.4b8t+00-3.3l6E-0  l"1  *3a7t+00- 
2. 5788-0 1-1 .3708-0 1-2. 5008-0 1-2.7578-02- 
2.710E-01  3.771E-0i.-2.S>J4E-01  7.512E-01- 
1  • 3908-0  l-2»4VSE+00~  1  •  3«;9E-0  1-1  » 4038+00 ' 
1  «0-i  lE~o  1-  1  .0078-0 1-i .OOOE-O 1-2.9992-02- 
l«084E~ol  3.8J0C.-0  a- i.  »  1  75E-0  1  7.6088-01" 
3.+77E-02-2.  sooE+00  3.323E-02-1  .4052+00 
2.576E-02-1 .bbbE-Ol  2.500E-02-2.739E-02 
2.7112-02  3. 0142-01  2.938E-02  7.604E-01 
1  .9878-0  1-2.492E-*  30  1  .6O0E-0  1-1  .3962+00 
1  .^eyt-ol-l .0418-01  1.250E-0 1-2. 31 38-02 
1 • 3E6E~o 1  3.7O7E-01  1.468E-01  7.5518-01 
3.570E-0 1-2.474E+00  2 • 986E-0 1-1 .386t+00 
2»32oE~o 1- 1 .478E-0  a  2 .^5oE-0 1“ 1 • 7558-02 
2.4398—0 1  3.6948-01  2.641E-01  7.459E-01 
4 • J59E— o 1 -2  #4622+00  3  .646E-0 1-1 .3782+00 
2.O35E-0 1—1 »442fc— Oi  2 . 7O0E-0 1 “ 1 .4442-02 
2.9608-01  3.6492-01  3.227E-01  7.400E-01 

5  .y  2bE-o  1  -2 . 4C9E+00  4  •  962E-0  1  - 1  •  3568+00 
3. B65E-0  1  - 1  •  36bt-0  1  3 . 7boE-0  1-6.0802-03 
4.061E-01  3.S47E-01  4.3ybE-0l  7.2b7E-01 
7.0898-0 l-2.ay7E+00  3 .y42E-0 1-1 .3372+00 
4 • 6 J7t— o 1  —  1 • 303E— 0 1  4.S00E-0 1-3.310E-03 
4.671 E— o 1  3.461E-01  5 • 267E-0 1  7.131E-01 
7.O56E-0 1 -2 . 373E+00  6 .592E-0 1 - 1 • 322E+00 
5 . 1 52E-0 1 - 1 • 260E-0 1  3 . OOOE-O 1 - 1 .300E-04 
5.409E-ol  3.4OOE-0i  3.B47E-01  7.037E-01 
1 »26o8+oO-2» 1608+00  9. J49E-0  1-1 • 1508+00 
5»t>bOE-ol-l  .0008-0^  btbbOE-Ol  1  *0008-05 
9.3498—0  1  1  •  1502+00  1.2608+00  2. 160E+00 
a.b02t-ol-3.220E-Ol  6.9b0E-0l-l .0008-05 
1.1958+00  6.694E-01  1.730E+00  1.6008+00 
9*bbOt-o  1-1  »000c.-03  9.550E-01  1.0008-05 
2.38OE+oO-7.6OOE-0i  1.6008+00-1 .0008-05 


4. 3958-0 1 -7 • 34 68-0  1-4 .061 E-0 1 -3. 658E -01 
■3.750E-01  2. 10 lE-02-3»B65E-0l  1.487E-01 
■4.962E-01  1 .3628+00-S.925E-01  2.429E+00 
•2  •  934E-0 1-7.51 28— 0 1-2. 7 10E-01 -3*77  IE— 0 1 
2.500E-01  2.757E-02-2.b78E-01  1.570E-01 
•3.316E-01  1 .367E+00-3.965E-01  2.468E+00 
■1 .  17bE-Ol-7»608E-ol-l»084E-Ol-3»830E-01 
•l«000t-01  2»9y9E-02-l «03lE-0l  1.607E-01 
•1  .3298-0  1  1  »403E+00-l  »b90E-01  2.4958+00 
2.9388-02-7.6048-01  2*71 lE-02-3*8 14E-0 1 
2.500E-02  2. 7398-02  2.578E-02  1.585E-01 
3.323E-02  1.405E+00  3.977E-02  2.500E+00 
1 .46bE-0 1-7.5512-01  1 . 356E-01-3.767E-01 
1.250E-01  2.313E-02  1.2b9E-0l  1.541E-01 
1.660E-01  1.3982+00  1.987E-01  2.492E+00 
2.64 l E-0 1-7. 4592-01  2.439E-01-3.694E-0 1 
2.250E-01  1.755E-02  2.320E-01  1.478E-01 
2.9bbE-0l  1.3862+00  3.570E-01  2.474E+00 
3.227E— o  1-7.400E-01  2.9e»oE-Ol  -3.649E-0  1 
2.7b0E“0 l  1.444E-02  2.83bE-0l  1.442E-01 
3.6462-01  1.3782+00  4.3b9E-0l  2.4622+00 
4 . 395E— 0 1 —7.2572—0 l  4.061E-01-3.547E-0 1 
3.750E-01  8. 0802-03  3.865E-01  1.365E-01 
4.962E-01  1 . 356E+00  5.925E-01  2.429E+00 
5.267E— 0 1-7.1 3  IE— 0 1  4.87 1 E-0 1 -3. 461 E-0 1 
4.500E-01  3.310E-03  4.637E-01  1.303E-01 
5.9422-01  1 . 337E+00  7.089E-01  2.397E+00 
5.8472-01-7.0372-01  5.409E-01-3.400E-0 1 
5. OOOE-O  1  1.300E-04  5.152E-01  1.260E-0I 
6.592E-01  1 . 322E+00  7.85BE-0»  2.373E+00 
7.4372-0 l-b.Sb3E— 01  6. 312E-01 -2.057E-0 1 
6.3122-01  2.0578-01  7.437E-01  5.553E-01 
l .7302+00-1 .800E+00  1  •  195E+00-8.694E-0 1 
6.9502—0 1  1.0008-05  8.802E-01  3.220E-01 
2. 140E+00— l .300E+00  1  .378E+00-4.643E-0 1 
1.3762+00  4.6432-01  2*1402+00  1.300E+00 
1.600E+00  1.000E-05  2.380E+00  7.8008-01 


2.500E+00-1 .COOE-03 


1 

4 

5 

10 

77 

88 

ay 

100 

173 

182 

183 

1  90 

166 

167 

177 

1  78 

34 

35 

46 

47 

130 

131 

142 

143 

2.500E+00 

l.OOOE 

-05 

1  1 

18 

19 

28 

101 

1  12 

113 

124 

191 

19o 

197 

200 

186 

187 

193 

194 

58 

59 

70 

71 

154 

155 

29 

40 

41 

52 

125 

136 

137 

148 

201 

202 

196 

199 

82 

83 

94 

95 

53 

64 

65 

76 

149 

160 

161 

172 

106 

107 

1  18 

119 

3*0808-01  3 . 0802—0 1  —  1 . 336E— 0 1  1.3368-01-7.2808-02  7.28o8-02-3.390E-02  3.3908-02 
3*fcOOE-o3-3.2oOE-03  3. 440E-02-3. 4408-02  5.0008-02-5.0008-02  6* 100E-02-6. 1 008-02 
6. 34 0E-Q2-6. 3402-02  6 . 360E-02-6»360E-02  6.3608-02-6.3608-02 
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UTRANL-II  Program 


With  the  present  mesh,  two  cases  of  flow  are  to  be  computed.  These 

are 

M  =  0.804,  k  =  0.253  (Subcritical) 

00 

M  =  0.903,  k  =  0.228  (Supercritical) 

00 


Input 


Listed  below  are  typical  input  cards  to  the  UTRANL-II  program  to 
compute  flow  field  for  the  above  cases. 

UNSTEADY  SOLUTION  —  NACA  64A006  —  202  NODES  —  M*0.B04  —  K*0.2S3 


.BOA  0«2bJ  1*1  0«25 

,o  o«  c« 


16*  32* 


0. 


(Insert  cards  generated  in  UTRANL-I) 

UNSTEADY  SOLUTION  —  NACA  6AA006  202  NODES 

•903  0*223  1«1  0*25  0«b 

•  0  0  •  0* 


—  K*0»903  K*0 • 223 

16.  32*  0. 


The  above  cards  indicate  there  is  no  non-zero  initial  guess  furnished  to  com 
pute  the  case  M  =  0.804,  while  the  case  with  higher  Mach  number  will  use 
results  computed  for  lower  Mach  number  as  initial  guess.  In  all  cases,  the 
nonlinear  boundary  condition  on  the  actual  airfoil  surface  is  to  be  used,  as 
suggested  by  IOPT(3)  =  0  in  the  options  card. 


Output 


The  output  from  this  program  includes  printouts  and  punch  cards  for 
the  solution  for  the  last  two  time  steps. 


The  printouts  are  in  the  following  order: 


•  Title  of  the  problem  under  study 

•  Convergent  limit  (a  positive  number  for  steady  flow) 

•  Control  keys  specified  in  the  options  card 

•  Reduced  frequency,  amplitude  for  the  oscillating  flap, 
x-coordinates  for  the  hinge  point  and  the  trailing  edge, 
free  stream  Mach  number,  initial  time  level  when  com¬ 
putation  is  started,  mean  angle  of  attack,  circulation  for 
the  airfoil  at  its  mean  position 

•  Total  number  of  elements,  number  of  nodal  points,  and 
the  full  bandwidth  of  the  resulting  matrix  equations 

•  Element  numbers  and  element  node  points 

•  Nodal  numbers  and  their  corresponding  x-  and 
y -coordinates 

■  Boundary  nodes  for  each  type  of  boundary 

■  Slope  for  nodes  on  the  airfoil  surface 

•  For  each  time  step,  the  computed  results  at  all  nodal 
points  are  printed,  together  with  a  printer  plot  dis¬ 
playing  the  pressure  coefficient  on  the  airfoil  surface. 

•  For  each  time  step,  the  unsteady  pressures  on  the  air¬ 
foil,  for  the  upper  and  lower  surfaces,  and  their  differences, 
are  printed.  At  the  end  of  each  period,  the  mean  pressures 
on  the  airfoil,  for  the  upper  and  lower  surfaces,  and  their 
differences,  are  also  printed. 


The  above  items  except  problem  parameters  and  computed  results  have 
also  been  printed  in  the  UTRANL-I  program,  but  are  repeated  here  for  com¬ 
pleteness.  For  brevity,  only  results  for  a  typical  time  step  are  listed  on  the 
next  five  pages  for  reference. 
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Printout  for  NACA  64  A006  atf  the  16th  Time  Step 

MA6H  NUMd£R=  .St*.  TIME  INCRtlENT  =  .7761  TIME  LEVEL  =  12.42 
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2.  FLOW  OVER  A  NACA  64  A410  AIRFOIL 


This  problem  was  analyzed  using  the  mesh  shown  in  Fig.  8,  which  con¬ 
sists  of  173  elements  with  202  nodes.  The  procedures  described  for  the  NACA 
64  A006  airfoil  problem  also  apply  here.  However,  the  data  coded  on  the  DATA 
statement  in  the  BLOCK  DATA  subroutine  should  correspond  to  NACA  64  A410. 
Also  modified  slopes,  taken  to  be  half  of  the  slope  at  0.5%  chord,  are  imposed 
at  the  two  leading  edge  nodes.  Input  cards  (page  66)  and  some  typical  output 
(pages  67  through  71)  for  UTRANL-II  are  listed  here  for  reference.  The  non¬ 
linear  boundary  condition  on  the  actual  airfoil  surface  is  used. 

The  case  to  be  computed  is: 

M  =  0.700,  a  =  2  deg  (barely  critical) 

00 

where  a  is  the  angle  of  attack. 
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Input  Cards  to  UTRANLi-II  for  NACA  64  A410 


i_1FTInO  SOLUTION  NACA  64A410 

1 

700  0.0  0.0  0.£5 

z  a.  o. 


20a  NODtS  —  M  =  0»  700 

0.5  1*  20. 


(Punched  cards  from  UTRANL-I). 


Results  for  the  last  time  step  are  listed  on  the  next  five  pages  for 
reference . 
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A  FORTRAN  listing  of  the  source  deck  for  the  UTRANL-I  program  is 
presented  in  pages  74  through  82.  The  program,  as  presently  configured, 
requires  46aK  words  to  execute. 


PROGRAM  MAIM  I  NPUT  =  65  .OUTPUT*  131  .PUNCH =65  *  TAPc.5*  1  NPUT  « 

1  TAP£.6=0uTPUT ) 

C  Gfc.Nc.RATc.  2-D  MESH  CONSISTING  OF  i  R I  ANGLES  ANO  UUAORlL.ATfc.RALS 

C  TO  INVOKE  ANY  OF  THE  FOLLOWING  OPTIONS  THE  CORRESPONDING  KEY 

c  I opt < i >  is  read  In  as  i 

C  I  OPT ( 9 )  — mesh  uENERATED  Is  FOR  AIRFOILS  OTHER  THAN  THE  0*06 

C  CIRCULAR-ARC  AIKFOlL 

C  I  OPT <  12) — PLOT  MtSH 

C 

DIMENSION  TITLE! 12)* I OPT (12) 

DIMENSION  NDELI400.4) .X<400)  «Y(400) 

DIMENSION  N1DSI 3) .N1D<3» 100) .VAF< 100) 

EuuIVALcNCE  INIDS< 1 ) » NFARF ) . <NIDS<2) .NWAKE ) «  <NIDS<3) .NBODY) 

NPM=400 

NEM=400 

NCN= 10 

C  READ  ANO  GENEKATE  DATA 

10  READ  <5.805)  (  T  I  TLE < I > , I  =  1 . 12 > 

IF  ( EOF  <  5 ) )  2000*20 
20  CONTINUE 

READ  820.  < I0PT( I ) . 1*1 . 12) 

CALL  MESHBW(NtM*NPM.NDEL.X.Y.NEL5.NP5.MAXN.N 'D.VAF.  I  OPT  ) 

READ  825.  NJDS 
DO  120  1*1 .2 
NS=NIDS< 1 ) 

120  READ  825.  < N I  D  (  I . J ) . U= 1 .NS ) 

NBw=6*<MAXN+l ) 

IF  <  I  OPT  <  I  2  )  *E0.  1)  CALL!'  MPLOT(NEM.NPM.NELS.NPS.NDEL.X.Y) 

C  PRINT  OUTPUT  DATA 

PRINT  910.  < TITLE! I ). 1=1.12) 

PRINT  820.  < IOPTI 1 ) . 1=1 ♦ 12) 

PRINT  920.  NELS  t  NPS . NBW 
PRINT  930 
DO  220  N= 1 .NELS 

220  PRINT  825.  N» ( NDEL < N . J ) . J* l . 4 ) 

PRINT  93 5 
DO  230  1=1  » NPS 

230  PRINT  940*  1.  X(l).Y(t) 

PRINT  951.  <NIDI 1. 1 ). 1*1 .NFARF ) 

PRINT  952.  (NIDI  2* 1  ) . 1* 1 .NWAKE > 

PRINT  953.  (NIDI  3. I ) . I*1.NB0DY) 

PRINT  955.  < VAFI I ) . 1* 1 .NBODY) 

C  PUNCH  OUTPUT  DATA 

PUNCH  825.  nEcS.NPS.NBW. (NIDS( I ) *  I* 1 .3) 

PUNCH  825.  (  (NDEL!  I . J) . J*1 .4) . 1  *  1 .NELS) 

PUNCH  840*  (XI  I )  »Y< 1  ) .1*1. NPS) 

DO  350  1=1.3 
NS=N1DS( I ) 
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JbO  PUNCH  82b «  (  N  ID(  I • U ) » U* 1 *NS) 

PUNCH  840 «  < VAF (  I ) *  1  *  1 .NBODY ) 

805  FORMAT  <12A6) 

820  FORMAT  (40 1?) 

825  FORMAT  (  1 6  I  «=5 1 
840  FORMAT  (1P8E10.3) 

910  FORMAT  < 1H 1 « 2A • 12A6// ) 

9*0  FORMAT  (1H0.-N0.  OF  ELEMENTS* 1  *  14 «  •  NO. 

I  i  FULL  tl  AND  WIDTH* 1  ♦  14//) 

900  FORMAT  ( ' QLLt *  NO.  AND  ELC.MENT  NODES'/) 


OF 


NUUE5= 


14. 


9J5  FORMAT  ('0  NODE*.  6X . ' X ( I ) ' . 1 2X. • Y (  I  ) ' / ) 

940  FORMAT  (16.  2E15.5) 

9b 1  FORMAT  ('ONODLS  AT  F APF I ELu ' / ( 20 l 5 > > 

9b2  FORMAT  ('ONODcS  ON  LINE  OF  SYMMETRY • / ( 20 1 5 )  ) 

*53  FORMAT  ('ONODt-S  ON  THE  A I RFO  I L  •  /  (  20  1 5  )  ) 

955  FORMAT  ('OSlOPE  ALONG  NODES  ON  A 1RF0 I L ' / ( BE  1 5 .5 ) ) 


GO  TO  10 
2000  STOP 


END 


75 


SUBRGUT  lNE  MESHOWl NEM  *NPM  ♦  NOEL.  •  X  »  Y  »NELS»  NP5  •  MMXN* 

1  N 1 O  * VAF  « !  JpT ) 

C  GENERATE  MESH  Ai\D  COMPUTE  THE  MAX t  D 1  FFt.Rt.NCt  IN  NOUta 

D 1  MENS  1  UN  NDEL(NtM*4>  *X(NPM)  *Y(NPM> 

DIMENSION  1 DUMP l 4 ) «  DUMP ( 4 ) 

DIMENSION  N I D ( 3  » 100) « VAF<  100 ) • i OPT ( IE) 

1020  FORMAT ( 1615) 

1 1 bC  FORMAT  (Jl5,5FlO»0*  15) 

c  generate  element  node  numbers  and  their  maximum  difference 

DO  100  N= 1 »nEM 
DO  100  U=  l  •  4 
loo  NDEL(N«J)=C 
NELS  =  0 
MAXN=0 
KUuNT=0 

IcO  Rc.au  1020.  NPt  #  N.E  I  *KE  J  •  KD  1  <  KDU  •  1  DUMP 
IF  ( NPE  .EO.  0)  GO  TO  1 29 
J1 =-KDJ 
DU  124  0=1 «KtU 
Ul  =  J1  KqJ 
J2  -  JI  -  KD I 
DO  124  1  =  1  • K  t 1 
NELs  =  NEuS  +  1 
Jc  =  J2  +  kd 1 
DO  124  K=1 « NPt 

124  ND£L(NELS*K)  =  IDUMP(K)  +  J2 

125  CONTINUE 

DU  126  1=1 • NPt 
DO  126  J=1 «NPc 

ND1FF=NDEE (NEoSt I ) -NOEL < NELS . J ) 

126  IF  (lABS(NDlFF)  »GT •  MAXN)  MAXN= 1 ABS1ND1FF ) 

GO  TO  120 

1 29  IF  (NELS  »Gt»  NtM)  GO  TO  300 
c  generate  nodal  coordinates 

READ  1020*  nPO 
READ  1200*  XLc  »XTE 
XCENT  = ( XLE+XT t )  *  *5 
FAcT=l./(XTE-ALt) 

1200  FORMAT  (2F10.6) 

KOuNT  =0 

l->0  READ  1150*  NO  »  NN  I  «ND  1  .RATIO*  DUMP*  1 D 
IF  (ND  *EU.  0)  GO  TO  140 
IF  ( 1O.EO.0)  GO  TO  132 
IF  ( KOUNT •NE»0)  GO  TO  131 

IF  (  I0PT(9> .Ew»l)  CALL  RAIF  ( I U • DUMP ( 1 ) « DUMP ( 2 ) « V AF S ) 

131  KDUNT=KUuNT+1 

NlQ<  3»KUUNT  )  =  .*D 

CAUL  GOCR  (10*  DUMP ( I )  • DUMP ( 2 ) *  DUMP ( 3 ) «  DUMP ( 4 ) « 

1  XCENT. FACT. VAFS) 

IF  (  1  OPT ( 9 ) • tU • 1 )  CALL  RAF  (  1 U • ( DuMP ( 1 )  +  •  5 )  *  I 00 • « DUMP ( 2 ) « VAF S ) 
VAF( KOUNT )=VAFS 
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132 


CONT 1NUE 
TEMP2  =  1*0 
IF  (NNI  tLT • 
jl  =  NNi  -  1 

TEMPI  =  1  «0 
TEMP  2  =  0«0 
00  134  I  =  I  • Jl 


3)  GO  TO  13b 


134 
1 3b 


Tt-Mp2  +  TEMPI 
kAT 1 0*T cMPl 

< JUMP l 3 ) -DUMP Cl))/ TEMP  2 
( DUMP ( 4 ) -DUMP ( 2 ) ) /TEMP2 


137 


140 


300 

olO 


400 

410 


+ 

+ 


TEMP3 

TEMP4 


TtMPi 

IEMP2 


TEMP2  = 

TcMP  1  * 

TtMP 1  = 

TEMP2  = 

Jl  =  ND 
TEMP3  = 

TEMP4  = 

DO  137 
X(J1)  = 

Y(  Jl  )  = 

Jl  =  Jl 
TEMP3  * 

TEMP 4  = 

TEMPI  * 

TEMP 2  * 

GO  TO 
IF  ( NPS 

WCTUPN  .  rv 

INSUFFICIENT  ARRAY  DIMENSION  assigned 

FORMAT3 MOTOo’mANY  ELEMENTS*  «bX« • NELS= • » 1 4  »  5X  » • NEM=  *  » I  4 

call  exit 

TJZSWo rooS;T  mooes-  sx-mps-.,*.**™  — 

CALL  exit 
END 


0*0 
0»0 

I  =  1 t  NN 1 

DUMP ( 1 ) 

DUMP! 2 ) 

+  no! 
a  TEMP3  + 
a  TEMP4  + 
a  WAT io* TtMPl 
a  WAT  1 0*  T  tMP2 
130 

•GT,  NPM)  GO 


TO  400 
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bUBPOUT  INC  MFl  JT  (  XC  «  YC  «  NODE  «  NN  •  NE  •  NELS  •  NPb  ) 

USE  CAlCOMP  p._orTc.K  TP  PLOT  GENERATED  McbM 

CODED  BY  CApT*  otPPY  VANKEUREN  OF  AFFDL.  WR  l  GHT-P  ATT  EPSON  AFb 

D  1  MENS  1  ON  XC  (  UN  )  »YC  (NN)  «  NODE  ( NE  «  4  )«XP(7)»YP(7) 

DEFINE  PLOT  SiZt  AND  PLACE  PEN  IN  INITIAL  POSITION 

DO  5  l=I»NPS 

XC(L)=XC<L>*8» 

S  YC(L) =YC(L)*8. 

c*m_l  plot < o.o »- j«o «-3 > 
call  Pl.uT  l  1  ,  0  »  3  •  0  •  ~3 ) 

DRAW  EACH  ELE^t^r  IN  TURN 
DO  10  I=1«NELo 
NPE  =  S 

IF  (  NODE  (  1  *4 )  • tO .0 )  NPE-4 

NPP=NPE-1 

DO  9  J=I  « NPp 

XP(  J) =XC(NOdE*  I  * J)  ) 

9  YP<  J)=YC(NODcU  »J)  ) 

XP(NPE) =XP(  I  ) 

YP(NPE) =YP(  I  ) 

XP( NPE+ 1 ) =YP( NPt+l  )=0.0 
XP(  NPE+<2  )  s't'P  (  NPt+2  )  =  1  »0 
10  CALL  LINE! XP« YPtNPt. I »0»0) 

white  nodal  nuMoEpo  and  finally  end  the  plot 

DO  dQ  L=l»NPb 
PL=FLOAT(L) 

EO  CALL  NUMOEPtXUL)  .YC(L)  ..  14.PL.0.-I  ) 

CALL  PlOTE 

return 

end 
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bUdWGUTiNt  GOCK  (!•  X.YO.XF.YF.XCW.FACT.UYUXI 
OuNcWATiNo  THL  tJODY  Gt.OMt.TWY  FOW  A  CIRCW-AK-AWw  A*kFOlL 

FY(X.H.H)=SQKT(K*K-X*X>-H 

DFy( X«Y«H>=-X/ (Y+H) 

DATA  R»H/4a  l8i6t>67«4a  1516667/ 

xc=xcw#fact 
W2=XC**2 
XO=X*FACT 
YO=FY(XO«W«H)*I 
DYDX*DFY(XO. YO*I »H)*I 
IF  (  ABSIDYDX)  ioLi  I  •E~10)  CO  TO  100 
IF  ( Q  YDXa  Gc • 0 • )  DYDX= 1 •£- 1 0 
IF  < DYDXaCTaO • )  GYOX*-  1  •  t-  10 
100  CONTINUE. 

yp=-i«/oyux 
YYX*Y0-YP*X0 
YYX2*YYX**2 
YPl =YP**2+ I • 

0ISC=YYX2*YP**2-YPI*< YYX2-W2) 

SD=SORT(DlbC) 

IF  (YP*I  al_T aOa  )  6U*-SD 
XF a ( SD-YYX*yP ) / y  P 1 
YFsS0WT(W2-xF*XF )*l 
Y0= YO/FACT 

xf=xf/fact 

yf=yf/fact 

IF  ( ABS( XO+a b ) aGTa I aE-5 )  CO  TO  120 
YO  =  0  a 
GO  TO  130 

120  IF  ( ABS( XO-a 5 ) aGTa l aE-5 >  CO  TO  130 
Y0=  0  a 

IUO  CONTINUt 
RETURN 
END 
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SUBROUTING.  RAiF  ( 10 ♦ X 1 ♦ YH • DYDX ) 

C  GENERATE  Y-COuKDlNATE  AND  SURFACE  SOuPE  FOR 

C  l>ATA  GIVcN  IN  ttLOCK  DATA  SUBROUTINE 

COMMON/AFG/XL  <  30 ) «  XU ( 30 ) • Y YL  <  30 )  *  Y  YU  <  30 ) *  NN 
DIMENSION  XX<  30) * YY<  30) *XP<  30) * YP<30> 

1 .YpD< 30*2) «YPL( 30) « YPU(30) ♦ YYO ( 30*2) 

2 « XXD  <  30 «  2 ) •  XPD<30»2) * XPL < 30 ) *XPU<30> 

EUUlVALcNCE  <  YP0( 1 ♦ 1 ) *YPL<  1 ) ) *  < YPO< 1 *2) *YPU< 1 ) ) 

1 ♦ ( Y YD (  1*1)  «YYL (  1 ) ) 

2*  <XXD<  1*1) «XL( 1 )  ) 

3*  <XPD<  1  ♦  1 ) *xPl ( 1 ) ) «  <  XPD<  1*2)  *XPu( 1 ) ) 

N2=NN+1 

CALL  GMPG  (NN*  XL*XPL *YYL* YPL ) 

CALL  GMPG  <NN»  AU*XPU*YYU* YPU) 

RETURN 
ENTRY  RAF 
DO  760  1=1 «N2 
1ARG= ( ID+3 )/2 
XX ( 1 ) =  XXD< 1 « I AKO ) 

XP( 1 ) =  XPD< 1 « I ARG) 

Y Y (  1  ) *YYD< I « I AKG ) 

760  YP< 1 ) =YPD( 1 « IARG) 

DO  800  U= 1 «NN 
UP  sJ+1 
XU  =  XP<  U) 

XJP=XP( UP) 

IF  ( (X1-XU)*(AI-XUP) .LE.O* )  DYDX  = YP ( U ) + ( X 1 -XU ) * ( YP ( UP ) - YP ( U ) 
1  /(XUP-XU) 

IF  (  U « Gc. •  NN )  oO  TO  8 00 
XXU=XX( U) 

XXUP=XX(UP) 

IF  ( (XI-XXU)*(XI-XXUP) •LE*0* )  YH  = YY ( U ) + ( X 1 -XXU ) 

1  *( YY<UP)-YY(0 ) ) /(XXUP-XXU) 

800  CONTINUE 
YH=YH**C1 

IF  (  X  I #OT • 1 • t —  3 )  RtTU«N 
YM=0« 

DYDX=0* 

RETURN 

END 
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SUBROUTINE  GMPG  (NN. XX « XP « Y Y * YP ) 

DIMENSION  XX(30> .yy«30).XP(30).yP(30) 

NlsNN-l 

N2sNN+  1 

DO  700  1= 1 »Nl 

IPs  1  +  1 

yp( iP)»tyy( iP»-YY(  i ) )/<xx« ip)-xx(  i  >  > 


XP( IP)*(XX( IP l+xxt  1  ) >  *.5 
700  CUNTINOC. 

XP( 1 ) *XX( 1 ) 

YP(  l ) =YP(2  > ♦( XP ( 1 ) -XP ( 2  > )*(YP 
XP(N2)*XX(Nn) 

YP{N2)*YPINN)t(XP(N2>-XP(NN)  ) 

RETURN 

END 


{3)-YP(2> ) / ( XP ( 3 ) —XP ( 2  )  ) 

* ( YP  C  NN ) -YP  C  N 1 )  ) / ( XP ( NN ) -XP l N 1 


)  > 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BLOCK  DATA 

DATA  FOP  NACA  64A006  AIRFOIL 

NN  NUMBeR  of  POINTS  DEFINING  THE  THICKNESS  DISTRIBUTION  ON 

EACH  SUkFACF  OF  THE  AIRFOIL.  _unM 

XL  ARRAY  SiJkInG  X-COORDINATfc.  IN  PERCENTAGE  «tASU«^  FHOM 

leading  euGe.  jf  nodes  on  the  lower  surface  or  tho  airfoil. 

XU  array  stoking  X-COOROINATo.  in  PERCENTAGE  MEASURED  FROM 

leading  edge*  of  nodes  on  the  upp&r  surface  of  thl  airfoil. 

YYL  array  STOKInG  Y-CuORDINATE,  in  PERCENTAGE.  OF  NODES  on 
THE  lOwoR  SuKFACt  OF  THE  AIRFOIL. 

VYU  array  stoking  y-coordinatl,  in  percentage.  OF  nodes  on 
the  uPPt-R  SuKFACt  of  the  airfoil. 

COMMON/AFG/XL i 30  » ’^O ( 30 ) . Y YL  t  30 ) . Y YU ( 30 ) . NN 

data  NN/26/  „  ,  „  J. 

data  (XL<  1  >, l *1.261/0... 5.  .7b.  I  .2o«2.5.5.  »7.t>*  ‘  ^ 

1  JO.. 3b.f 40. t AS. tbO.«bb..tjO. ,ob..70..7b..b0..BJ..y0*.^^*» 100»/ 

data  <xu<  1  )  ,l«l*2«»)/0...t>.  .75.  1.2b.2.b.b.*7.a«l0.  .  lb.  .20.  ,*.b.. 

1  j0..3b..40..4b.*b0..b5..t>0..bb..70..7b..B0..ob..90..GB..100./ 

DATA  ( YYL( I ) . I  =  lf20)/0. .-.465. -  .  bb5  ,  -  .  739  ,  -  1 . 0 1 6  ,  -  1 • 399 . 

1  -l.6a4.-1.9lY,-2.cbJ,-2.bb7,-E.7S7,-2.c^6.-c.977.-E./99. 

J  -2,S(4B.-E.a23»-2.Pb3.-2.4  3b.--E.  I  db  *  -  k  .907.-1  .602,-1 .26j. 

3  -.967«-.64gf-. 331 .-.013/  . 

DATA  < YYU( 1 >. 1=1. 261/0... 4B5..S6b,. 739, 1.016. 1.399. 1.664. 1.919. 

1  c.  2bJ,  e.  bb7. 2. /  =  7, 2. 696.2.977,2.9*9, 2. *4  =  ,  .  o2b ,  6.  bb  J  .  2 . 4  Ja  . 

2  2.  166.1.907.1.602.1. 265. .967 , .649, . 33 1 , .013/ 

END 
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o  r 


dLOCK.  DATA 

sTnSsjj**  toT,\i  ;ssiu  x*  on 

^CH»£S;*S5o!IUM5-Ji22!«Tt*  1N  ""“TTS- 

XU4Ui“A^STj.U^NX-cL«UlN^^0lN^^^XAGLr«tASUH^  F*D» 

™  ss:  ssas  of  *»„  «* 

.^3oU“xXU,30 , 

‘  i;:^.«:to.W.^6V.«v«:.7..tt7*.7V.0.V.S*.MA.0V.b'«. 

is:si;!^;;jU;!J4ts:foS.TO.u*.«...b..-b..A..«>.»<- 

. . 

rH:oS::t:n:::i::?s:-t:o::::t^i:;s:-::;t::^:-:sn: 

«j  o.624.o.l06.a.^.4.7o.J.v67.J.0U»w:.0JO,l.0<:o..04:l/ 

End 
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uuu  O  U  U 


c 

c 

c 

c 

c 

c 

c 

c 


MAIN  1|,.POt.6».i~*l»UT*»S.-HUNtH.6S..I~<A3.lNMUT. 

-----  :--A^^T^«eNTS 
OF  PHEVIOUH  CA*E  AE  ATARI  IN6 
•JHILE  THt  OTHtR  OPTION  IS  IGNORED 

1  OPT  (  4  )  =  1  »  UPDATING  THt  CIRCULATION 


COMMON  /COM . / AM (3 . '.Hulas, . *. zoa » . 
nc^;^t^:3%,;vAFK,,o, 

DIMENSION  TlTc£(12>.IOPT>l£).NIDSO>.AH<eO> 

ECU  l  VAUtNCC**?  MI  l)S<  1  >  «NFAHF  >*(NIDS(2>  *NWAFE  >  •  IN.  US'  3.  * NODDY  } 
L)mTm  gamma/ 1  *40/ 

NPM=202 

NEM=173 

NRM=NPM*3 

NCMs84 

NdCT*50 

EXP=-GAMMA/ ( GAMMA- 1 «0> 

CONST*0«5* ( GAMMA- 1 .0) 


read  title* 


CONTROL  keys*  and  problem  parameters 


100  READ  (5*845)  ( T  I  TLE< I  >  » I X1 » 12) 

IF  ( EOF ( 5 ) )  2000  »  105 
105  CONTINUE 

HEAD  *  83E^°hmAC^KAEC  »  ANOC*  MNP  •  XTE  •  DTM  *  TP  TO  *  aTEeT  *  ALPHA  «C.  SCO 


READ  835 t 

TITLE,  CONTROL  KEYS,  and  PROBLEM  PARAMtTtRb 


PRINT 


PRINT  910*  ( TITLE! I) ,1*1  * 12> 

PRINT  B22*  i<SC»ANDG»XHP»XTE*RMAC*TF IN»T0» ALPHA. CIRC0.2TEST 

1 RES* - 1 
IPRO*DTM 
SUMAC«RMAC**2 
RKSC«RKSC*2. 

ANdGs  AN0G*3«  1 4  i  5926/  1 80  * 

ALPHA® ALPHA#3 • 1415926/ 180* 

IF  (RKSC.GT.O*)  0TM=6.283I853/<DTM*RKSw) 

NlFTsTFIN 
T 1 = ( TO— 3* )*OTM 
T2®  T 1 +DTM 
T3s  T  2+DTM 

IF  (  I  OPT ( 1 )  *t0*  1>  GO  TO  382 


C 

c 

c 


READ  AND  GENERATE  MESH  INFORMATION 


READ  825*  NELS.NPS*NBW*<nIDS  1  »I«I»3> 

REAO  825*  ( ( NODI l*0)*J®l»4)»I*l *NELS ) 

READ  840*  lxm*V«*»»l*l,NPS* 

DO  110  1 = 1 «3 
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_  ■  Cm 


o  n 


N3  =  N IDS  <  I  ) 

1»C  PlAd  b2b.  <NlU( 1 , J) ,J=l,NS) 
~--A0  b4C.  (  vAi-  (  1  )  ,  1  =  i  ,NuODY) 


C 

C 

C 

c 


IF  ICPTIj 
AND  apply 

ImONl  1  NEAa 


)  =  1*  Mb Yd  NODES  ON  A  1 KK OIL 
lI.mlAnI/liD  bOUNuAPY  CONDI 
DOUaDmPY  CONDITION  oN  Tnt 


SUPKACc  TO  CHGWDLINt  POSITION 
TlON.  OTrtLK^lSL.  11/POSl 
AIWFGIl.  SUWFaCl 


IF  (  I  OPT  (  3  )  .  iME  •  1)  GO  TO  lib 

Nu.VtN  =  NoooY/2  *2 
Ju  1  IS  J=  1  ,NcV'cN»2 
*  =N ID  <  J  .0 ) 

V(  I  ) =-0.015 
I  =  N1D<  3.J  + 1 > 
lib  Y< 1 )=  0.015 
1  I  6  LJiNir  1  Not 


paint  GoNokaTco  jata 

NFt9w  =  Nfcw/b 
NtQ=3*NPS 

cOO  PAINT  920.  NtLb.NPS.NbW 
PP I  NT  9 JO 
po  220  im=  1  .  NtL  J 

220  PPJNT  625.  N* ( NOD < N.J J ,J= i ,4 ) 

PHI NT  93b 
DO  230  1=1 .NPS 

tOC  POINT  940,  1,  Xli).Y(l) 

PP 1 NT  9b 1 «  ( ft ID(  1 ,  1  )  ,  1= l ,NFARF ) 

Paint  9b2,  <  ,\i I d<  2 » l » »  1  *  1  .nwakc ) 

PAINT  933,  3, 1  )  ,  1= 1 .NbUOY) 

PH 1NT  9bb,  ( VAF (  1 )  «  1  a  i  , NoODY ) 

DC  244  11=1 ,NPS 

P'V'L  (11)  =Pi»i  Ac 
244  CUNT  1 Not 

c  A-AJ  IWLHO  INITIAL  LLtSS  W  PHULEtU  » 1 T«  &H!)  aol_UT10N 

IF  (  1  OPT  <  2  )  .ImE.  1)  GO  TO  250 
auAD  84C,  < SLP (  I  )  ,  1  a  1  . NLO ) 

Pi- AD  04  u  ,  <  S  I  I  .  NJA/  )  ,  1  =  i  ,  NtU  ) 

Co  TO  2*3 

llC  DO  260  l  =  l»Ntvj 

-jt  1  (  1  )  =U. 

toe  SLP ( 1 ) =0.0 


c. 

c 

V- 

c 


INTlGPATIoN  bTAKT  Ptwt  and  CheCK  IF 
PP1NT  Fail  to  Co.WlPGC  ANu  PPOLotD  TO 
C ON T 1 NUt  TO  I  NT  CoWATE 


N  1  FT 
NlXT 


IS  LXCctUcD.  if 
CASL.  OTFitPw  I  Sl 


265  T1=T2 
T  2  =  T  3 
T o  =  T3  +D  T  >1 
IkES=1PcS+I 

F  oRiNuLATc.  oYoTtN  OF  ALGtoPA  l  C  cQUATICNS 
DO  266  1=1, NCO 
OD  c66  0=  1  ,  Nd.v 
tbfa  D (  1 «  0 ) =L. 0 


SO, 
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(..All  I  o^G  (  0  •  *  bw  vIal  ) 

(,«Ll  <  ->G,«1hC  »  p<RN»NCM»NLG*Nu*  •  New  •  NEI >  »  NcU  »  aLP  •  3  »  COF  tNWvi*X  ,  Y  • 

1  LL  1  ) 

(.  IMPOSE  G.G,  CgR  r«KhlL.LUi  LlNL  Of  5  YMMl  T  R  Y  «  ANu  ON  aIrFCIL 

L«LL  L'UY_ON  (  I  Ui  »  N  IL)  »NoLT  «VAF  «  X»  Y  »NPN.  »  5  »  aLP«  NFo'i,  NCM, 

1  Nglit  Ml)*  *  AL|  HA  »  k/iAL  «  L  I R CO « KEY  •  XHP,  ANUo  *  aX3>C  «  i>L  I  *  I  OPT  > 

L  LALL  UANOeO  bdLVLH  TO  SOLVE  THe  SYSTEM  OF  EGUATIONS 

Call-  dNOc.0.  ( 3  «  l»X. 'i  «NCP,  «  NEo  *NHb* ) 

C 

0  Pk  I  NT  Co y,PUT c. o  RESULTS 

cVb  PRINT  9  70 *  PmAC,0TM»T3 

b  ALL-  TLyA\/(l#*JO'*lAC) 

CALL  OJTP  ( 6 « SEP, 6L I »RML* COF  •  SuMAC » OgNST *cXP» 

I  mNOG  «  I  Ana  «N3P  ♦  ■'iO'H  *  Nw.v, «PCTE»  Nl  D  *NdOl)Y  «NbCT  «  I  PRO ) 

PR I  NT  9ob 
1PNT=C 
I 6TCP=0 

NeVEN  =  NcJuOY/2*2 

ue  -320  g = l « nE V EN ,2 

IF  <  J,Ge,Neven~1  )  ISTOPsJ 

I =  NlU< 3*0) 

IF  ( I0PT<4).t-.i)  GO  TO  3  15 
bP=S<  J*3) 

CALo  F  PLlT  ( 60  »  I PNT  »  AFtLN » 1 STuP  »I»I»X(I)*CP) 

Go  TO  320 
.5 lb  CONTINUE 

CP  =  S  <  I  «  I  ) 

call  FPLUT  < bo »  i PNT ,AR,LR»  ISTOP. I .2. XI  I  )  .CP) 

I  =NI0<O» 0+1 ) 

CP  =  S  < I ♦ 1 > 

L/VlL  FPi_oT  (  60  »  1  PNT  ,  AR»Lk  «  I  STUP  *  2  *  2  »  X  (  I  )*  CP  ) 

J20  CONTINUE 

c 

o  CHECK.  CUNVeRoeNce.  IF  SO,  PUNCH  THE  FINAl  SOLUTION  ANp 

C  PkOCEEO  To  NeXT  CASE •  OTHERWISE,  UPOATe  SOLUTION  ANO  CONTINUE 

l.  TO  INTEGRATE 

C 

IF  ( IRES. LT.l. Ok. PCTE.GE«ZTEST>  oO  TO  362 
IF  <  IRES.oE.N 1FT  )  GO  TO  600 
GO  TO  700 
6C0  PRINT  9«0 

700  PUNCH  040,  ( SLp (  1  )  *  I  =  I »NEG ) 

punch  640»  < s( I .row) » 1=1 »neu> 

Ir  (  I  OPT ( 4 ) „ eg • 1  )  PMINT  325,  CISCO 
Go  TO  100 
362  COnT  I  NOe 

UO  390  I = I »Nto 
SH  (  I  )  =3LP<  I  ) 

390  b(_P<  I  )  =3  <  I  *NBa  ) 

C 

C  IF  I  OP  T ( 4 )  =  1  ,  UPJATE  The  circulation 

c 

IF  <  I  OP  T<4)«t.L*«0)  GO  TO  265 
IL=NID<2» 1 )*i-2 
lU  =  NI0<  2,2 )*3~2 
lIRlO  =  SLP<  I u ) ~ S  LP  <  IL) 

PRINT  325,  CinCO 
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C  1  klUlA  T  1  ON  =  •  *  c.  1  4  •  b/  ) 


GO  TO  2 6d 

j co  format  ( iho.-^omPuTed 

bdO  FORMAT  (40  12)  _  w  .  w  ,  r  . 

0£;2  FORMAT  (lrto,-  *n5C  =  '  .F9.5.4X.  'ANDO  =  •  *F9.a.  =X.  XrtP  =  ,F^.j. 

1  Jx  *  '  XTc.  =  •  «  F  9  •  3  »  3X •  •  WMAC  =  •  .  F  V  •  5 . 3X  .  •  TF  1  N  =».F9.b.3X* 

2  »T0  =*»F9.5/1H  »'  ALPHA  =  •  « F  9.  b  •  3X  «  •  C  I RCO  =  '»F9.5*3X* 

3  • ZTEST  = • »  F9 • 5/ ) 

025  FORMAT  (  16  l  5) 

( 8F  1  o  •  0  ) 

( 1P8E10.3) 

(  1H  ,1 3A6 ) 

( 1  rl  l  •  1/A6/  » 


835 

840 

b4S 

910 

9l0 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


NO*  OF  N0JcS='«l4» 


(IrtO.-NO.  OF  ELEMENTS=* • 14. • 

*  POLL  bANO* IOTM= • « 14// ) 

('OElE.  nJ*  and  element  NODES'/) 

(.OolD  NODE  •  «  10X.'X(1)'.12X.'Y(1)'/) 

(  I  10 . 2E 15.5) 

(•ONCDtS  AT  FARF 1 ELO • / ( 20 15)) 

(•ONoDES  ON  LINE  OF  SYMMETRY •/( 2015) ) 

(•ONOOES  ON  TPtE  A  IRFO 1 E ' / ( 20  1 5 ) ) 

COSLOPE  ALONG  NODES  ON  AIRF01L'/{8EI5.5)  ) 

( 1 H 1 i  -MACH  NUMbER= ' »F6.3* 

1  5X  .  '  T  I  Me  INCKuMENT  =  '  «  F6*4  »  5X  «  •  T  1  ME  i_EVEE  ='«F6.2/) 

S»dO  FORMAT  COFAlL  TO  CONVERGE  IN  SPEClFlcD  NO.  OF  ITERATIONS') 
905  FORMAT  (  1H  1  ) 

2000  STOP 

end 


930 

935 

940 

951 

952 

953 
935 
970 


FORMAT 

format 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


SUBROUTINE  OUTP  l  S  ♦  SLP*  SL  1  «RME  «  CGF  «  SUMAC.  «  CONST  «  EXP  « 

1  AnUG«  I  RtS  »NPS  #  RBW  *  NRM  «  PCTE  ♦  N  IU  tNbOOY  »NdOT  »  I  PRu  ) 

PRINT  CUT  RESULTS  FOR  THE  CuRRc-NT  TIME  STC.P 

UUMMON  /CuM  I  /  mM  I  3)  •  AMT  (  3)  •  AMTT  l3)*Tl«T<>«T3»UT 
O  I  MENS  I  UN  S  t  ,\KM  ♦  RBW  )  ♦  SLP  I  NWM  )  •  ^L  1  I  NWM  )  iRMl.  (  NPS  )  «COF  <  NPS  ) 
It  Nl  U  (  3  tNtiCT  ) 

PRINT  y75«  IKES 

IND-O 

PCTE=C* 

DU  305  U= I t NPs 

I  I s3*J 

I  ids  1  1  —2 

POT  =  S  t I EtNbw ) 

UPT  -  S  <  I  1-1  tNOtf ) 

V  =  S<  I  I  *  NOW ) 

U-UPT  + 1 «0 
USQ=U*U+V*V 

PHT*POT*AMTi3)  +  SLPt  I  2 ) * AMT ( 2  >  +SL 1  I  I  2  >  *AMT (  1 ) 

ASQ-CONST*  (  1  •0-USU-2t*PHT  )  ♦ 1  • /SUMAC 
RMxSORT(USU/ASU) 

PRA= ( 1 .+CUNsT*RM**2 ) **EXP 

CPs*-2«*UPT 

CPU=-2**PHT 

CP=CPS+OPU 

CPU = CP 

IF  ( AfcJSt ANUG) tGTt 1 »E-b)  CPU  =CP/ANDG 
S( J  « 1 )=CPU 

CGF<  J)=SUMAc*l  1  •  0+2 •  4*UPT  > 

UELMsRM— RMc ( J ) 

PADWs Aos I DELM ) *  100  t 
RMl<  J ) skM 

IF  ( PAOMt LE # P CT E )  GO  TO  305 

IND*U 

pcte*paom 


CP. DEEM 


JOb  PR l NT  97b.  J .PO T »U • V .COF ( J ) •  RME ( O ) < PWA .CPU. 

PR I  NT  d  1  0  *  InD.PCTe 
INpRU=  IRtS  / I  PHD# I  PRD- 1  RES 

PRINT  905 
NEvEN=NOODY/2*2 
DO  500  I  =  1  .NcVEN.2 
11=1+1 
1 E  =  N I U ( J •  1  ) 

I U  =  N I U ( 3  «  1  1  ) 

DCP=S(  I L  » 1  J -3 (  I U  » I  1 
S< I .3 J=DCP 

IF  ( INPRD.Ne.  0>  W  TO  400 
S ( 1 « 2  J  =(S(I»2)  +S(  lE.l J#.5)/IPRD 
3(I1.2)=(S(ll»2>+S(!U.l)*.5>/I  PHD 
D5uM=S(  I *2 J-a (  I  1 #2) 

PR I  NT  908.  S<  Il»2> «3(  !*2J*D3UM 
b( I .2)  =S(  IL» 1  )*.5 

S(  I  1  .2 J  =S(  Iu»l )*.& 

GO  TO  500 

400  a ( 1*2)  =a(I,2)  +3(lL.lJ 
S (  I  1 «  2 ) =S(  I  1  • 2 )  +S  ( IU. 1 ) 

500  PRINT  990  .  I  *  IU.  IC«S(  lU.  1  )  »S(  IE,  1  J  .DCP 

BIO  FORMAT  (  1  HO  <  —NODES  =  •  .  13.5X.  ‘MAx(UEcMJ  PCT  *  •  .  IPi_  10.3/  ) 

975  FORMAT  ( 1H0.4X.-T1ME  STEP  ='.14/ 

1  lHO  . 3X  » 'NODE- »7X» 'PHIT • « 10X. »UCOM»  . 10X.  'VCOM* ♦ 1 IX.  'COF • . 10X. 

2  'LMAC . 1 1X.-PRA-.1 IX. ‘CPU ' . 12X. »CP» . 10X. • DEEM • / ) 

97b  FORMAT  ( I O. I P  9h  a  +  «4 ) 

9o5  FORMAT  (  1H0. 3X  I  •  .faX*  •  Io»  .OX*  •  IC»  .  7X  .  •  CP-UPPCH •  .  7X  . 

1  'CP-L04eW‘ , 7X.-0CP*  »  9X .  'AVE-CP+* .QX. • AVE-CP-* .aX. • STEADY -DCP •/ J 
9bb  FORMAT  (75X.1P3E15.4J 
990  FORMAT  (31 l0.1Picl5.4) 

RETURN 

end 


90 


-  jw&m  -  .•(iKnir  I  1-1  - 


-jU(3i«0U  T  iNt  Bl)Y  CUM  MUb.NIU«NdCT  «  VAF  .  X  .  V  .  NPM  .  b  ♦  5LP  «  N*M  .  NCM  . 
1  NLU  «  NOW  ,  ALPHA  *  KMAL  «  C  lWCOiKLYi  XHP  •  AN  Jo  •  Wf,jC  •  OL  1 , lUPT ) 

IMPOSITIONS  JP  THt  BOUNDARY  CONDITIONS 

COMMoN  ZoOMi/hMI  3)  •  AMT  (  3  )  •  AMT  T  (  J  )  «  1  1  *T2*TJ*dT 
DIMENSION  NIUM  3  1  «NID<  3«NBCT  )  ♦  VAf-  (NUCT  >  ♦ X < NPM >  «  Y  (  NPM  )  « 

1  S  (  NkM  «  NCM  )  ,  S  oP  (  \AM  )  oL  1  <  NWM  )  *  I  OPT  (12) 

DATA  P 1/3.141^926/ 

NO* l =  NBW- 1 
NMQW=NdW/2 
NFAkf-  =NlOS<  1  ) 

NPA IR=NlDS<  2 ) / 2 
NdOD Y  =  N I Dd ( 3 ) 

CNsT=0.b*ClRCO/PI 
NOOT=SCNT( 1.0-RMAC*WMAC) 

DO  120  0=1 «NoOOY 
I  L.3-N  I  □  (  3«  J  ) 

It  =  3* 1 1 J 

OUAN^VAF  (  o  )  -ALPriA 
XN£DU=X< 1E3)-XHP 
DDUAN=0« 

IF  ( XREDu*CT»0» )  00  TO  90 
ARGU=KKSC*T3+1. 5707963 

DouAN  =  An0G*( cost APGU>-XREDU*  RKSC*S1N( ARGu > > 

IF  (  XRE.OU«CE«  1  .t-5>  DCiuAN=DGuAN*.b 

90  continul 

IF  ( 1CPT(3> .tu.l)  GO  TO  96 
DO  95  K=2«N0*/i 

95  M lt-1 «K) *S( It-i#K)+OUAN*S( IE»*-1 > 

S(  I E- 1  «N3W ) =S (  I t-1 «NB*/>  +UUAN*S(  It  *NbW> 

96  CONTINUE 

DO  1 00  <= 1  « Ndw 
100  S( IE *K) =0*0 

S  <  I E  «  NHBW )  =  1 , 0 
S( It .  NBW) =QuAM+00UAN 
IF  ( 10PT<3> ,tU.l)  GO  TO  120 
S ( I E . NHUW- 1 ) s-OUAN 
lc:0  COnT  INUt 

CALC  tdwv  ( I • ) 

DO  220  N=1«nPmIk 

1c=J*NIU(2«2*.i-1)-3 

lUs3*NlU<2«r*W)-3 

NC*IC-Iu+NHLW 

DO  210  M=2«3 

IE: 1 U+M 

I tL  = 1L+M 

ID= Iu-1L 

10P= ID+1 

DO  190  K*  1  DP  tiYbMl 
190  S  (  I  EL  «  K ) *S (  I  El  *  X, )  +S  (  It*  K—  I D  > 

S( ItL«NBW)=S( i£L»NBW)+S( lt*Nb*> 

DO  200  K= 1 «Nbft 
dOO  S(  1 E  «K  > =0 • 0 

S(  It*  NHOd/ )  =  1  »0 
G(  IE*  NC ) =  —  1 »  0 
IF  (  M •  NE •  2  )  GO  TO  210 
S< It.NHBW-1 ) =AMT ( 3) 


S< 1E.NC-1 ) =  -AMT ( 3) 

S(  I E ♦ NbW ) =  AmT ( 2 ) *( SLP( I EL— 1 ) 
1  -bi-1  < lu-l  >  > 
clo  CONTINUE. 

220  CONTINUE 

DO  340  N= 1 «nFARF 


-SUP (  1  £-  1  >  > ♦AMT ( 1  >  * ( SL I 


l=NlU(  l *  N ) 


J=3* 1-3 

YM00  =  ^0OT* Y  <  I  ) 

WSQ=YMOD**2+X l I ) **2 
YY= -YMOD 

THETA  =  ATAN2 < YY *X  l  I  )  ) 

IF  (THETA  «lT*  0.0)  THET A=THETA+2. *P I 

DO  310  K=1  *  3 


IE= J+K. 


DO  300  U=  1  .f> 

300  S( IE.U) =0*0 
310  S< I E • NHBW ) = 1 *3 

S(J+1»NBjO)=  ChST*THETA 
3<  j  +  2.N6\m ) =cNoT* YMGD/PSU 
S<J  +  3.NfcU))=-OST*P00T*X<  I)/«SU 
340  CONTINUE 
4lTU«N 

end 


1  EL - 1  > 


one 


\ 


SUflWOUT  I  Nt  FPDO  T  (  Ml  f  1  PNT  *  AW  «I_W  •  1  STOP  i  NC  .NCMAX  «  V  1  «  V 2  ) 

WeSULTS  PoOtT  i N<J 

LUGICAL  i-M(  2)  *LN(  120)  .eP«4  )  *L.X(4>  *t-R«  1  )  •  PC  <201 
lNTLGt.w  nI(2)«ST<2)«SG(2) 

D  I  Mb.  N  b  I  UN  AR(M1  )  «N(  120  >  «WHO(  30  )  *Z(  lb)  •  l  OFF  ( 2  >  « 1 5P <  2  > • SF 1 ( 2  >  » 6F2 <  2 ) 
j  ,  I D  (  2  )  *  1  sC  *  2  )  *1  1  ( 4  )  •  1P14)«1M(4)*IX(4>.  IP  <202  >  •  1  A<  23b> 
tuul  VAi-t-NCS  (  IP  (  1  )  «UP(  |  )  )  •  ‘  IM(  1  >  «l_M(  1  )  )  *  (l_N(  1  >  .N(  |  )  •WHO!  1  )  )  • 

I  (  I  X  (  I  )  .LX(  I  )  )  •  (M2C.W1  (  I  >  »  •  (N58.RI  <2>  ) 

JATA  (Z(l)«l  =  l*i!>)/l..l.<:b«I.3.1.75«2».2.5.3.«J.3»4.«4.5.3.«6.« 

1 7. , d. .9*/ 

2 AT  A  (|sP(I)«1=1«2)«(IX(I).I=1«4>«<1P(I>*I»1*4>*<IM(I>.I=1.4)/ 

I  lC.5«4*lHX.4*lH+»4*0/ 

OAT A  10 ( I ) . I Dl 2 ) « IBG( 1 ) . 1SG( 2)/ I  20 .38* I •  -!/ 

OAT  A  LN‘  I  )  M  1  »LCU  )  *LC<2>/1 •bU*5W.O. 1  UHL*  1  HO/ 

IK (  I  STOP) 173*172.172 

172  J- 1 PNT  +2 
I  F ( J.GT.M1  )G0  TO  173 
IPnTs J 
An  (  J-  1  )  *V  l 
AW ( O ) *V2 
DO  10  K*  1  *  4 

10  1M(K)  *  NC 

J  =  J/2 

LW(  J)  *LM(2  ) 

173  IF(NC.LT.NCmAX|Rc.TUHN 
IK ( 1 STUP.oO. 0 IKtTUWN 

C  1  0  3  We AO ( s« 1 )LNi 1  )  *  (  I  I  <  l  )  «  I  *  l  «  4)  «  *t-C<  I  >  ♦  l*  1  »20)  •  <  LN<  1  )  •  I  *2*52  > 

1  FORMmT ( A1 t 2  I  3 »2 l 1*71 Al ) 

4W|T£.(6*  1  >LN<  1 )  •  <  I  I  «  I  >  ♦  1*1  *4)  ♦  < i-C C  I  )  ♦  1*1  *20)  ♦  <LN<  I  >  .  1*2*52  > 

OU  171  1=1*2 

W I (  I  >  * ID(  1  > 

SG(  I  )  =  I SG(  1  ) 

1F(  I  1 (  I ).Nt.0)Rl  (  1  )  =  I  l<  1  * 

IF ( I  I ( I +2 ) «EQ. 1 ) SG( l > =  -SG< I ) 

171  ST(  I )  =  («!  <  I  )  +  l-SG(  l )*<RI < I >-l) >/2 
1F( WI ( 1 ) .GT. 120 ) R 1 < I) =1 20 
1F( Nl (2) .GT .2 J8) *1 <2) =238 
N20  =  N 1 20/6 
N 1 20  =  6*  N20 
DO  3  I  a l • 2 
AMI N  =  AW< I ) 

AMAX*AR< I ) 

DO  7  J*I » IPnT  *2 
A  A’.  =AH(  J) 

IF  1  AM IN.GT «  aA  1  ) AM  IN* A A  I 
7  IF < AMAX.oT . aAI ) AMAX=AA1 

IF  f  ABS ( AM  1 N )  .LT*  *000001)  AMIN  =  0*0 
IF  X  ABS  <  AM  Ax )  •  e  T  •  .000001)  AMAX  =  0.0 
W  =  AM'AX-AMI  N 

IF  (  AbSlW  )  «lT  «  1  •  t-9»  AND.  AtoSI  AMAX)  *OT*  1  »t-9  )  R=1.0 
IF  ( AdS(W) .LT.l.t-y. AND. AMAX. LT.O* )  W*-AMAX 

IF  ( AdSI W ) *lT • 1 . t~9 *  AND* AMAX.GT .0*  >  W*  AMAX 

DO  22  J=1  «  15 

B=AL.UG10lW/(rtI  1  I  )-2)/2<J)  ) 

M=0 

IF ( L)  .LT .0 ) M=M— 1 


I  1 

12 

13 

14 
13 
16 
17 


19 

20 
21 
22 

23 

24 

25 

26 
27 
2b 

29 

30 

31 

32 

33 


36 

37 
3d 

39 

40 

41 

42 


43 


47 

48 

49 

50 
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' 


; 


_ —  ■ 


c  =  z< J)*10<*» ( M+l  )  51 

L  =  AM  IN/ 2  (  J  )/  1  J  .  »  »  C  M-t-  I  J  c,2 

11=0  53 

IK ( d.LT.O) 1 1 = 1 1- 1  54 

IK  (  (fill  I  1-2+1  1  |  IN)  1  tt  <  1  v  »  1  V  5b 

IB  0=10. *C  56 

19  IK  (  J.tu.  1  )  SM  1  iuI.  57 

22  1K< 0.LT.3M1M oMiN=C  SB 

5F] (  1  )  =  < l./SMlN)*bG(  1  )  59 

b=A^lN/BMlN  60 

M  =  B  61 

IF ( B.LT.O >MsM-  1  62 

5F;><  1 )=ST<  1  )-M*SG<  1  )  63 

KfHOC  I  )  =  SF 2  (  1  )  +0  •  5  64 

M=SK2(1)  65 

UU  25  J=1  •  10  66 

IK  (  M-J-l  (M-J  )  /  I  5P(  1  )  )  *  I  SP(  1)  ) 25.3*  25  67 

25  GONTlNufc.  66 

3  1  OFF (  I 1=0  69 

UU  1 C 1  1=1 <Nbd  70 

10  1  1  A(  1  )  =0  7I 

DO  102  J= 1 , IFNI ,2  7 

1 w ( J ) =SF 1 ( 1 ) *AR ( J ) +RHO ( 1 )  73 

1 T  =  5F 1(2) * AR I J+ 1 ) +RHO  C  2 )  74 

1F( J.NE. 1 ) GO  TO  109  75 

lk< 1PNT+2) -Z  76 

1W<  J+l ) =0  77 

106  1 3= 1 T  7B 

GO  TO  102  79 

109  I F ( IT-131104. 103»105  felO 

104  1K( J-M  )*ik( 1PNT+2)  a  1 

1W( IPNT+2)=J+1  Uc 

GO  TC  106  b3 

10b  1  =  1  T  +  1  a* 

106  1=1-1  as 

1 1  =  1  A ( 1  )  56 

1F<1D1O6<1O6<107  87 

107  ih<  J+i >  =  i«< 1 1 )  ae 

iw< l i >=j4i  e9 

102  1  A  <  1T)=J+1  90 

LAST  *  1 PNT  +  2  91 

J J=  1  OFF ( 2 )  92 

lZh=5F2<I>  93 

L2V=SF2<2)  94 

DO  100  l*l<Nbd  95 

DO  40  J= 1 <  NJ  20 
40  N< J)=  1 H 

1F< <  1-1 )*<  1-N58 I  I 140< 151 < 140  96 

151  DO  141  J=  1  <  Nl  20 
14  1  N< j)  =  I X  <  1  ) 

140  LN<  1  I =LX(  1  )  10  1 

lMN1201=LX(  1  )  102 

1F<  LZH.l_t.0.0K.L2H.GT<N120  >  GO  TO  131  103 

LN(  LZH)  =L.P<  4  )  104 

1 J I  NB= 1  1 05 

1F< 1 .NE.JJ1G0  TO  35  106 

N5=2  107 


I  3=  I  OFF (  1  )  10* 

UJ  32  J=13.n120«10  11C 

3d  UN ( J ) =  LP(  4 1  111 

IF(  1 .NE.L2V1G0  I J  35  112 

00  13b  J=1  «Nl20 
13b  Nl J 1  =  IP(  1  ) 

35  13= 1 Al 1 )  11= 

I F  (  I  3 .  EU •  0  )  GO  T U  1  2  1  Ho 

ICQ  LAST= IW(LAaT)  1 17 

12  = [RluAbT-i )  llO 

1  1  *L  AST / i  11® 

UM(  d ) =LH(  1  1  )  120 

IMM  =  1M(2) 

LN ( 1 2 ) =LC ( INM)  121 

IF< LAST.NE. I  3 ) GO  TO  120  122 

121  GO  TO  (38.41 ) . Nb  123 

38  wHJ TE< 6.3* 1 < N< J ) .3* 1 «N120> 

3*  FORMATl 1 lrt  • 120A1 ) 

GO  TO  100  126 

41  AA 1*1  127 

VAluE*( AAl -SF2 < 2 ) ) ^6F 1 < 2)  12b 

WHITE (6.42 1 vAuUfc . (N<  J) . J* 1 .N120) 

42  FuRMATllH  • 1 P L 1 0 • 3 « 1 20A 1 ) 

ICO  CONTI  Nut  131 

13= 1  OFF (  1 1  132 

J=0  133 

DO  49  1  =  13  «Nl 20  *  10  134 

J^J+l  1 3b 


AA=  l 

49  RHOl 3)  =  ( AA-SF2 ( 1 ) 1/SF1 ( 1 ) 
1 F (  1  OFF ( 1  1-5)62.62.63 

62  WHiTE(6.50)  (  RllO (  I  1  .  1  =  1  •  J  1 

50  FORMAT!  10X. 121 1PE10. 3)  1 

HE turn 

63  1F( J.GE. 12 1 J=  1 1 

WHITE (6.64 1 ( RHO ( I ) . I = 1 . J 1 

64  FORMAT (  16X. 1 1 ( IP 1 10.3  1  1 
RETURN 

END 


EgjHSHWEBIM  .Jim 


oUBRUU  T  l  NE  Nt  * K.  (  SUMAC  .  NWM  *  NCM  *  NE U  .  NBW  •  NE M  •  NtLi.  NvjL)  •  SOP  *  s  • 

1  COL>  •  NPM  *  X  »  Y  ♦  SC  1  > 

c  GENERATE  SYSTu-M  MATRIX 

DIMENSION  COL>-  (NPM)  «X(NPM)  .Y(NPM>  *  XU ( 4 >  *YU(4> 
COMMuN/COMS/Aut  12*  12)  .6U( 12. 12) .AS(9»9>  *BS(9.9> 
COMMON/COM9/CS  JA  (  12  >  *  CSOB  (12)  *CSLA  (  9  >  •  CSLB  (  9  > 

D  1  MENS  I  ON  NOD  (  NEM  «  4  )  •  S  (  NRM  «NCM  )  ♦  SCP  (  NWM  )  ♦  BE  1  (  NRM  )  •  PMD  (  12*2) 
NHBW 

JO  4BC  N= I «NcLS 
11=1 

IF  ( NOD ( N • 4 ) )  402*402*404 

40  2  NP£L  =  3 
NTRS=I 
GO  TO  410 
404  NPEC=4 
NTRS=4 
12  =  0 

DO  408  1=1*4 
N I  =  NOD ( N  •  1  ) 

406  IF  (COEF(NI)  »GT •  l .00)  12=12+1 

IF  (12  *£0«  0)  NTRS=2 
IF  ( N •  CE • 5 )  NTRS=4 
4  10  DO  425  1  =  1  •  NPtL 
N I  =  NOD ( N  » 1  ) 

XUl  1  )=X(N1  ) 

YU< I )=Y(Nl ) 

DO  425  0=1 *3 
1S=J*(NI-1 )+J 
1E=J*( 1-1 >+J 
PMD( IE* 1 )=SLl( I 0) 

425  PMD( 1E«2)=SlP( Is) 

CAlc  EMoT  (  XO»  YO  «PMD«  SUMAC  •  NTRi> ) 

00  450  1=1 1 .NPEL 
NM=  3* ( NOD ( N , I  ) -1  ) 

1 E=  3* (  I-l ) 

DO  450  11=1.3 
NR=NR+ 1 
1 E= 1 E+ 1 

S(  NR  .  NBA  )  =  6  (  NR  *  NtiW)  +CSOA  (  IE  )  +CSQB(  1 E ) 

DO  450  0= 1 .NPtL 

NC=3+ (NOD( N, J ) - 1 ) -NR+NHBW 

JE=J*( J-i ) 

DO  450  00= 1 , 3 

NC=NC+1 

0E=0E+1 

450  S ( NR.NC ) =S ( NR  »  NC ) +AU (  1E.OE)+BO(  It.OE) 

460  CONTINUE 
RETURN 
END 
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5UBW0UT  1NE  EMUT  t  Xu*  YU  *PMUD  «  SUMAC  *NTM5 ) 
oENtWATC  MATKiX  FOW  A  UUADR  1  LATERAL  OH  TW 1  A  NOLL 
CUMMON/CoMb/ Au  (  12*  12)  *BU  (  12*  12)  *  Ao  (  S' *9 )  toi(V«V) 
CJMi«ION/LjM<v/CiJA  (  12  )  *  CSOt!  (12)  •  CELA  (  y  )  *  CBLB  (  9  > 

L)  1  MLNj  I  ON  XQ(  4)  «  YU ( 4 )  *  XT  <  3  >  *  YT  (3)  *  MP  (3*4) 

O  1  Mt.No  1  ON  PMULH  12*2  )  *PMTD<  9  *  2  ) 

DATA  MP/  1*2,J«3»4*1*2*J*4*4*1*2/ 

FTOW- 1 .  J 

1 F  (NTNS  • LQ •  4)  FT0P=0*5 

CO  ICO  1=1 « 12 
CSG A (  1  ) =0* 

CSQS<  l ) =0  « 

DO  100  J- 1 * 1 2 
AU<  1 . J)=0* 

BU<  1  . J ) =0* 

1O0  continue 

NTRA  = 1 

IF  <NTRS.EU.1.QR.NTPS.EU*4)  GO  TO  210 
NTRA  =  2 
NTRS=4 
DO  200  1=1*4 

IF  (YU(l).CT.  l.E-5)  GO  TO  200 
NTrA  = 1 
NTRS=2 
GO  TO  210 
200  CONTINUE 
210  CONTI NUt 

DO  150  1  lsNTRAn'lTWS 
DO  105  1=1*3 
N 1 = MP  <  1 . 1  I  ) 
l Ts  3* (  1-1 ) 

10  =  3* (Nl-l  ) 

DO  102  0=1*3 
1T= 1T+1 
10= IU+1 
OO  1 02  <= 1 *2 

102  PMTD(  |T*<S)  =FMuD(  1U*X) 

XT  <  1  ) =XU(N1 > 

105  YT  <  1  ) =  YU  <  N 1  ) 

CAUL  EMTC  (XT *YT*PMTD* SUMAC) 

DO  130  K=1 *3 
NW=3*(MP(K« I l )-l> 

1£=3* ( K- 1 ) 

DO  1 30  KK= 1 . 3 
N«=NR+ 1 
1E= IE+1 

CSQA ( NW )  =  CSQA ( NR ) +CSLA ( IE)*FTOR 
CSQB ( NR ) =C  SQb ( NK ) +  CSLB ( 1E)*FT0R 
DO  1 30  C= 1  «  3 
NC=3*<  MP(L  « 1  I ) - 1 ) 

UE=3*(L-1 ) 

DO  130  LL= 1 ,3 
NC=NC+ 1 
J£= JE+1 

AU<  NR.NC  )  =  AQ(  iMR  *NC  )+AS(  I  E  *  OE  )  *FTOR 
BU(NR,NC)  =BG( NR*NC)+BS<  IE*  JE)*FTOH 
130  CONTINUE 
150  CONTINUE 
return 


SUBROUTINE  FMTC  < XL. YL.PtLU. SUMAC) 

EVALUATt-  tLciMCNT  MATRIX  FOR  A  TRIANOLE  BY  OAUSS1AN  UUADRATuWE 
COMMON  /COMj /AM( 3 )  *  AMT ( 3)  . AMTT<  3) . T 1  *  T2»  T3.DT 
COMMON/CUMS/Aj (  12* 12)  «  DO <  1 2 »  1 2 >  ,  AS ( 9 ♦ 9 )  « BS C 9 ♦ 9 ) 
CUMMON/COM9/CSUA ( 12) .CSOBI 12) »CSLA<9> «CSLB<9) 

DIMENSION  DnX (9)  *P<  9)  «Q<9)  «NP<5>  »B<  3>  «C<  3)  *XL(3) . YL  C  3)  .S(3)« 

1  GpA ( 9 )  *DNXX ( 9  >  * DNY Y  <  9 ) . SN ( 9 ) »PELD<9»2) 

DIMENSION  ElM  (  J  ,7)  .  WT(  7) 

DIMENSION  S ?( 2) «W2(2)  ,PD<2) .UD<2> 

DATA  MMAX/2/* o2/-*577Jb027. . 57735027/ » *2/2* * 5/ 

L  AT  A  LMAX/7/  »<VT/0*22S,3wO*  1  32394  lb*  3*0  •  1  25939  1  8/ 

CATA  E I  NT/ j*0 • 3 A 333333. 0*0596 1 587 . 3*0* 4  70 1 4206 *0*05961587. 

1  3*0*47014205 .0*0596 15B7.0*  79742699 « 3*0* 10128651  . 

2  0. 79  74269 y*3*C. 10128651 .0*7974  2699/ 

DATA  nP/ 1 .2. 3 *1 *2/ • GAMMA/ 1 .40/ 

DC  2  1=1,9 
CSlA( 1 ) =0. 

CSLB ( 1 ) =0* 

DC  2  J= 1*9 
AS ( I « J) =0* 

5S(  I »  3 ) =0 • 

2  CONTINUc. 

DO  4  1=1,3 

U  =  NP< 1+1 ) 

N=NP< 1+2) 

to( I )=YL< J)-YL(K) 

4  C( I ) =XL ( K ) -XL ( J ) 

AWEA=0*5*<B<2)*C(3)-B<3)*C<2) ) 

VOLUME=APtA*DT 
CST 1=1 .C-SQMAC 
CST  2  =  SUMAC* ( 1 • O+GAMMA ) 

OO  100  L=1 .LMAX 
DO  10  1=1.3 
10  SCI) =E I  NT (  I ,L > 

CALL  DEWV ( AREA  *  d  »  C  «  S , DNX , DNXX  »  ONYY  »  SN ) 

0  =  0  • 

UX=0. 

DO  JO  1=1.9 
PLLD12=PELD( 1,2) 

U  =  U+DNX<  I  ) *PElD I  2 
30  UX=UX+DNXX(  I  >*PeLD12 
ALPHA=CST1-CST2*U 
DO  40  1=1.9 

Pll) =ALPHA*qNXX (  I ) +DNYY (  I  ) 

40  L<  1  ) =P( 1  )-CsT2*UX*DNX< I  ) 

DO  45  1=1.2 
PD< I ) =0* 

UD ( 1 ) =0* 

DO  45  J= 1*9 
PLLDJ1=PELD(  J»I  ) 

PD( I ) =PD ( 1 ) +PtLD J I *SN ( J ) 

UD<  I  )=UD(  l  ) +PELO J l *DNX ( J  > 
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4S  CUNT  1  Nut. 

DO  100  M  = 1  »  mMAX 

call  TDRV  <  S2 l M )  *  SQMAC ) 

WE.IG  =  WT(L)*W2  (M)*VOLUME 
DO  48  J= 1 ♦ 9 

48  GPA  (  J  )  xP<  J  )  -  (  AMT  t  3  >  *DNX  (  J )  -f  AMTT  <3>*6N< J)  ) 

SUM  =  0 • 

DO  50  J=l»2 

bO  6UM=6UM  +  ( AMT  (  J ) *UD  (  J ) ♦AMT  T ( J ) ( J ) ) 

DO  60  I  =  1 »  9 

CSA=-WEIG*  ( AMT ( 3  >  *ONX (  I  ) +AMTT ( 3 ) *6N ( I ) ) 

C8B=WEIG  *  Q|I| 

CSLA<  I  >  =CSLA (  i  ) +CSA*SUM 
CSLB< 1 )=CSLa( 1 )+CSB*SUM 
DO  60  J= 1  *  9 

AS  (  I  i J) =  AS (  I  • J ) +CSA*GPA<  J) 
uS(  1  * J) =38 <  I  ♦ J  J+CSB*GPA( J) 

60  CONTINUE 
100  continue 
return 
end 


SUBROUTINE  TDNU  (QPN. SQMAC) 

COMMON  /COM i /AMI  3) • AMT < 3 > . AMTT < 3 >  .  Tl  .T2»T3»DT 
RMA=S0MAC*2. 

RM3=RMA 
DT 1 =T2-T 1 
DT2=  T  3-T  2. 

DT=DTl+OT2 
DNl*l./IDT*r,Tl) 

DN2=-l./(DTi*DT2) 

DN3=1 •  /  <  DT *dT  2 ) 

RETURN 
ENTRY  TDRV 
T=( T 1+T3+DT*UPN ) *«5 
AM( l )«CT-T2)*( T-T3)*DN1 
AM( 2 ) * ( T“T 1 ) *  I T-T3)*DN2 
AM<3>=<T-T1 )*{ T-T2)*DN3 
AMT( 1 )=<T*2.-T2-T3)*DN1*RMA 
AMt(2)=<T*2.-T3-T1 >*DN2*RMA 
AMt( 3)  =  <  T*2.-Tl-T2>*DN3*RMA 
AMTT< l >=DN1*RMB 
AMTT<2)=DN2*RMB 
AMTT<  3) =DN3*RMB 
RETURN 
END 

i 
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SUBKOUT 1 Nt  DbKV ( AWL A  tb  « C  *  S  «DNX  * DNXX .DNYY . SN ) 
fc.VAL.UATt  THE  DERIVATIVES  OF  SHAPE  FUNCTIONS  AT 
DIMENSION  B  (  3  1 *C(3>  »S(3> «DNX  <  9  >  .DNXX ( 9 ) • ONYY ( 9 
DATA  Np/ 1 *2. 3  *  I ♦ l/ 

TvdQA-E»*AWEA 
T*OAbU=TWOA#*2 
H  =  S< 1 >*S<2)*S( 3) 

DO  200  1=10 


J  =  NP< 1+1 > 

K=nP< 1+2) 

S  1  3  S  (  1  > 

SJ=S<  J) 

SKaS(K) 

01 =B( 1 ) 

0J  =  B<  J ) 

0Ks  B ( K ) 

C1*C(  1  ) 

CJ=C ( J) 

CK=C(K) 

S1SU=SI*S1 
fcl  1 SQ=B 1 *b  1 
ClSU*CI*Cl 
ALFA*0«5*(CK-CJ) 

BETAaO«5* ( B  J— cJK. ) 

HX=B1*SJ*SK+BJ*S<*SI+BK*S1*SJ 

HXX  =  2 . * ( S 1 *B J*BA+S J*BK*B I +SK*B I *BU ) 

HYY  =  2.*<S1*CJ*CK.+SJ*CK*C1+SK*C1*CJ) 

CSS*6.*S1*« 1 «-Sl  J 
CS  =  6 • *  1 1 .-2»*SI ) 


L  =  3  *  1  “  2 

SN<  l.  )  =  Sl SU*  <  3  »~2  «*S 1  >+2.*M 
ONX( L  >  =  B  1 *CsS+2 • *HX 
DNXX <  L  >  =B1 SO*tS+2.*HXX 
QNYY  <  L>  =C1 so*cs  +  2.*hyy 
cs=ck.*sd~cu*si^ 


A  GAUSSIAN 
«NP( 5) »SN< 


L=L+1 

BN(  L  )  =SliU*Ci>  +  ALFA*H 

ONX ( L ) =2 • *B 1 *5 I *CS+TWOA*S 1 SQ+ALF A*HX 
DNXX ( L ) =2. *B I SQ*CS+4 . *B I *TWOA*S 1 +ALFA+HXX 
DNYY  <  L ) =2. *C I SO*CS+ALFA*HYY 
cS=BJ*S<<'—ts|^*SD 


200 


300 


-=L+1 

iN<  L) =SlSO*BS+BtTA*H 

DNx (L)=2»*B1*SI *BS+BETA*HX 

DNXX  <  L ) =2«  *B I SG*ttS+BETA*HXX 

DNYY  <  L ) =2 . *C I SU*bS+4 . *C 1 *TWOA*S 1 +BETA*HYY 


DO  300  1=1«9 

DNx (  1 >  =DNX<  I  ) /TtfOA 

DNXX ( I ) =DNXX < 1 ) / TWOASG 

DNYY(  I  )=DNYY(  I ) / TWOASQ 

WETUWN 


END 
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SUBROUTINE.  BNOEUI  A. NRMAX .NCMAX .N •  ITERM) 

EQUATION  SOLVtK  FOR  BANDED  NON-SYMME TR I C  SYSTtM  OF  tUUATlUNS 
USING  GAUSSIAN  ELIMINATION 

OULU T  JON  sTcRtD  IN  The  COLUMN  A<N»2*lTEKM> 

DIMENSION  At NRMAX .NCMAX) 

CERO  = I .E-06 
pare  =  CLRO**2 
NBND=2* I  TERM 
NoM  =  NciND  -  1 

BEGINS  ELIMINATION  of  the  lower  left 
DO  IOOO  1=1.  N 

IF  (  ABS( A ( I » I  TERM) )  •  LT •  CERO)  GO  TO  410 

GO  TO  430 

410  IF  (  ABS ( A (  I  » 1 T wRM ) )  .LT •  PARE)  GO  TO  1600 
PRINT  420.  A  (  I  t  1  TERM )  .  1 

4^0  FORMAT  <34H  WARNING — ILL— CONO I T 1 ONEO  A-MATR I  X . E 1 6. 6  * 5X *  I 4 > 
430  IP! T* I+ITEWM-1 

JLAST  =  M l NO (  IP1T.N) 

L  =  ITERM  +  1 
DO  500  J  = I .  JLAST 
L  =  L  -  I 

IF  (  ABS<A(J,L>)  .LT •  PARE)  GO  TO  500 
B  =  A< J.L) 

DO  450  <  =  L.  NriND 
450  A(J.K)  =  Alj.K.)  /  B 

IF  (I  .EG.  N)  GO  TO  1200 
50C  CONTINUE 
L  =  0 

JFIRST  =1+1 

IF  (JLAST  ,LE.  1)  GO  TO  1000 
DO  900  J=  JFlKSTt  JLAST 
L  =  L+  I 

I TL= I TERM-L 

IF  (  ABS ( A ( J  » I TL ) )  .LT*  PARE)  GO  TO  900 
JMl  =  J— L 

DO  600  KsITeRM.  N8M 
«ml  =  k.-l 

600  A<  J.KML)  =  A(  jML»KJ-A<  J.KML) 

A ( J  «  NBNO )  =  A  <  JML  »  N6ND ) - A ( J . NBNO ) 

NlT*N— ITERM+1 

IF  (I  .GE.  NIT)  GO  TO  900 
DO  800  K= 1 .  L 
NBNDK=NBND-K 

600  A  <  J  «  NBNOK ) s  — A ( J  * NBNDK ) 


900  COnT I NOt 

1000  continue 

c  bACK-SUBST 1TUTI0N 

1  tlOO  L  =  I  TERM  “  1 
OO  1500  1=2,  ft 
DO  1500  J=l,  C 
1 T  J= I TERM+J 
NP 1 =N+ 1 - 1 
NPl J=NPl+J 

IF  (NPIJ  • Gt •  N)  00  TO  1500 

A1NPI iNUNO) =A( NP I »NQND>-A(NPI J ♦ NbND ) *A < NP l ♦ ITU) 
1500  CONTINUE 
RETURN 

C  PRINT  THE  ENTIRt  MATRIX  IF  ZERO  ON  MAIN  DIAGONAL 

1600  PRINT  1601 

1601  FORMAT  <22H  IcRJ  ON  MAIN  DIAGONAL) 

DO  1602  1*1,  ft 

1602  PRINT  1603,  (A(I,J),  J=  1  »  NBNO ) 

1603  FORMAT  < 10E1 2.4) 

STOP 

END 
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